This file and all other referenced in the code can be found at the repository: https://github.com/Andros-Spica/Quaternary_Angourakis-et-al-2021

Load packages:

require(reshape2)
## Loading required package: reshape2
require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.6.3
require(patchwork)
## Loading required package: patchwork
## Warning: package 'patchwork' was built under R version 3.6.3
require(colorspace)
## Loading required package: colorspace
## Warning: package 'colorspace' was built under R version 3.6.3
require(grid)
## Loading required package: grid
require(gridExtra)
## Loading required package: gridExtra
## Warning: package 'gridExtra' was built under R version 3.6.3
require(png)
## Loading required package: png

Define lookup functions to access data between tables using identifiers or ‘keys’:

getRelatedVariableWithKey <- function(keyInCallingTable, keyInCalledTable, variableInCalledTable)
{
  return(
    sapply(keyInCallingTable, function(x) variableInCalledTable[match(x, keyInCalledTable)])
  )
}

getRelatedVariableWithTwoKeys <- function(firstKeyInCallingTable, secondKeyInCallingTable, 
                                       firstKeyInCalledTable, secondKeyInCalledTable, variableInCalledTable)
{
  temp <- 1:length(firstKeyInCallingTable)
  
  for (i in 1:length(temp))
  {
    temp[i] <- variableInCalledTable[firstKeyInCalledTable == firstKeyInCallingTable[i] &
                                       secondKeyInCalledTable == secondKeyInCallingTable[i]]
  }
  
  return(temp)
}

Define the format of image files to be generated (only “png” and “eps”):

listOfImageFormats <- c("eps", "png")

Crop table

Load crop table:

cropTable <- read.csv("models/cropsTable.csv", skip = 4)

Clean empty columns (generated with automatic names starting with X):

cropTable <- cropTable[,names(cropTable)[grep("^(?!X).*$", names(cropTable), perl = TRUE)]]

Reorder crops for a more logical display:

WARNING: doing this manually makes the script sensitive to the value of crop-selection in NetLogo. Thus, this chunk must be updated if crop-selection is changed before running experiments.

cropTable$crop <- factor(cropTable$crop, levels = c("proso millet", "pearl millet", "rice", "barley", "wheat 1", "wheat 2"))

Get vector of colours to represent crops:

cropColours <- qualitative_hcl(nlevels(cropTable$crop))
# alternative using only base R graphics
#cropColours <- rainbow(nlevels(cropTable$crop), s = 0.8, v = 0.8, end = 0.85)
knitr::kable(cropTable)
crop species cultivar T_sum HI I_50A I_50B T_base T_opt RUE I_50maxH I_50maxW T_heat T_ext S_CO2 S_water sugSowingDay sugHarvestingDay root.zone.depth..mm.
wheat 1 Triticum aestivum Yecora Rojo 2200 0.36 480.0 200.00 0 15 1.24 100 25 34 45 0.08 0.40 330 105 1000
wheat 2 Triticum aestivum Batten 2150 0.34 280.0 50.00 0 15 1.24 100 25 34 45 0.08 0.40 330 105 1000
rice Oryza sativa IR72 2300 0.47 850.0 200.00 9 26 1.24 100 10 34 50 0.08 1.00 170 320 400
barley Hordeum vulgare hypothetical; based on DSSAT v4.7 1762 0.42 350.0 170.00 0 15 1.24 100 20 34 45 0.08 0.40 330 100 1000
pearl millet Pennisetum glaucum hypothetical; based on DSSAT v4.7 - MLCER047.CUL - Middle Variety 1220 0.25 245.0 120.00 10 33 1.90 100 5 35 47 0.08 0.05 200 320 1000
proso millet Panicum miliaceum L. hypothetical 1328 0.29 157.3 96.75 0 18 3.00 100 130 34 45 0.08 0.05 107 220 1000

Experiment 1

Experiment 1 simulates the growth and yield of crops included in ‘cropTable’ using the NetLogo implementation of SIMPLE crop model (Zhao et al. 2019) adapted as a module of the Indus Village model under the default parameter configuration, which is aimed to approximate the conditions in Haryana, NW India.

Loading and preparation

Load data:

yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")

Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):

yieldData$yield[is.na(yieldData$meanARID_grow)] <- NA

Force the same order of crops used by cropTable:

yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))

Calculate the range of ARID and productivity (yield):

minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)

minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)

Load daily simulated weather and ARID data:

weatherData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_weather_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")

As the reference real-world data, we use the data downloaded at NASA´s POWER access viewer (power.larc.nasa.gov/data-access-viewer/) selecting the user community ‘Agroclimatology’ between 01/01/1984 and 31/12/2007. The exact location used as a reference is:

  • Rakhigarhi, Haryana, India (Latitude: 29.1687, Longitude: 76.0687)

We selected the ICASA Format’s parameters:

  • Precipitation (PRECTOT)
  • Wind speed at 2m (WS2M)
  • Relative Humidity at 2m (RH2M)
  • Dew/frost point at 2m (T2MDEW)
  • Maximum temperature at 2m (T2M_MAX)
  • Minimum temperature at 2m (T2M_MIN)
  • All sky insolation incident on a horizontal surface (ALLSKY_SFC_SW_DWN)
  • Temperature at 2m (T2M)

Load Soil Water (or Soil Water Balance) model code to calculate ARID given real-world data:

source("FigWeather_input/soilWaterModel/estimateETr.R")
source("FigWeather_input/soilWaterModel/watbal.model.R")

Load and format the real-world dataset so that they can be used in the Soil Water model:

weatherData_rakhigarhi <- watbal.weather.file("FigWeather_input/POWER_SinglePoint_Daily_19840101_20071231_029d17N_076d70E_5b401917.csv", 
                                              year = 1984:2007)

Estimate reference evapotranspiration using FAO recomendations (Penman-Monteith method):

weatherData_rakhigarhi$ETr <- estimateETr(weatherData_rakhigarhi$I, 
                                          weatherData_rakhigarhi$T2M, 
                                          weatherData_rakhigarhi$Tmax, 
                                          weatherData_rakhigarhi$Tmin, 
                                          weatherData_rakhigarhi$T2MDEW, 
                                          weatherData_rakhigarhi$WS2M)

Calculate ARID in real-world data:

Use Soil Water model to calculate soil water content proportion or WATp (mm mm-1 day-1) and ARID coefficient.

The first step is to set all soil parameters, here assumed to be constant (the same values used in experiment 1 simulations):

watbalParameters <- watbal.define.param()

Run the model inputing parameters and dataset:

watbalOutput_rakhigarhi <- watbal.model(watbalParameters, 
                                        weatherData_rakhigarhi,
                                        typeOfParValue = "nominal")

watbal.modeliterates for every day in the dataset calculating soil water content, absolute and proportion (WAT, WATp), and the respective ARID coefficient.

Show parameter values used here in Soil Water model (both in watbal.model and in the SIMPLE crop model). In experiemnt 1, they should be the same:

knitr::kable(
  cbind(
    parameters = c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF"),
    SIMPLE = as.vector(unlist(yieldData[1, c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF")])), 
    watbal = c(200, # these values of elevation and albedo are the default in estimateETr
               0.23, 
               as.vector(watbalParameters[1, c("CN", "DC", "z", "FC", "WHC", "WP", "MUF")]))
  )
)
parameters SIMPLE watbal
elevation 200 200
albedo 0.23 0.23
CN 65 65
DC 0.55 0.55
z 400 400
FC 0.21 0.21
WHC 0.15 0.15
WP 0.06 0.06
MUF 0.096 0.096

Summaries per day of year

Create data frames containing the daily summary statistics of simulations and real-world data:

  1. Create names vector (used for both datasets)
weatherSummaryNames <- c("dayOfYear",
                         "solarRadiation.mean", 
                         "solarRadiation.sd", 
                         "solarRadiation.max", 
                         "solarRadiation.min", 
                         "solarRadiation.error",
                         "temperature.mean", 
                         "temperature.sd", 
                         "temperature.max", 
                         "temperature.min",
                         "temperature.error",
                         "maxTemperature.mean", 
                         "maxTemperature.max", 
                         "maxTemperature.min",
                         "maxTemperature.error",
                         "minTemperature.mean", 
                         "minTemperature.max", 
                         "minTemperature.min",
                         "minTemperature.error",
                         "temperature.lowerDeviation", 
                         "temperature.lowerDeviation.error", 
                         "temperature.upperDeviation",
                         "temperature.upperDeviation.error",
                         "precipitation.mean", 
                         "precipitation.max", 
                         "precipitation.min",
                         "precipitation.error",
                         "ARID.mean",
                         "ARID.max",
                         "ARID.min",
                         "ARID.error")
  1. Simulation data
weatherSummary <- vector("list", length(weatherSummaryNames))
names(weatherSummary) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential

for (day in 1:365)
{
  weatherSummary$dayOfYear <- c(weatherSummary$dayOfYear, day)
  
  tempData_thisDay <- weatherData[weatherData$currentDayOfYear == day,]
  
  # solar radiation
  weatherSummary$solarRadiation.mean <- c(
    weatherSummary$solarRadiation.mean, 
    mean(tempData_thisDay$solarRadiation, na.rm = TRUE))
  weatherSummary$solarRadiation.sd <- c(
    weatherSummary$solarRadiation.sd, 
    sd(tempData_thisDay$solarRadiation, na.rm = TRUE))
  weatherSummary$solarRadiation.max <- c(
    weatherSummary$solarRadiation.max, 
    max(tempData_thisDay$solarRadiation, na.rm = TRUE))
  weatherSummary$solarRadiation.min <- c(
    weatherSummary$solarRadiation.min, 
    min(tempData_thisDay$solarRadiation, na.rm = TRUE))
  weatherSummary$solarRadiation.error <- c(
    weatherSummary$solarRadiation.error,
    qt(0.975, 
       length(tempData_thisDay$solarRadiation) - 1) * 
      sd(tempData_thisDay$solarRadiation, na.rm = TRUE) /
      sqrt(length(tempData_thisDay$solarRadiation)))
  
  # temperature
  
  ## daily mean
  weatherSummary$temperature.mean <- c(
    weatherSummary$temperature.mean, 
    mean(tempData_thisDay$temperature, na.rm = TRUE))
  weatherSummary$temperature.sd <- c(
    weatherSummary$temperature.sd, 
    sd(tempData_thisDay$temperature, na.rm = TRUE))
  weatherSummary$temperature.max <- c(
    weatherSummary$temperature.max, 
    max(tempData_thisDay$temperature, na.rm = TRUE))
  weatherSummary$temperature.min <- c(
    weatherSummary$temperature.min, 
    min(tempData_thisDay$temperature, na.rm = TRUE))
  weatherSummary$temperature.error <- c(
    weatherSummary$temperature.error,
    qt(0.975, 
       length(tempData_thisDay$temperature) - 1) * 
      sd(tempData_thisDay$temperature, na.rm = TRUE) /
      sqrt(length(tempData_thisDay$temperature)))
  
  ## daily max
  weatherSummary$maxTemperature.mean <- c(
    weatherSummary$maxTemperature.mean, 
    mean(tempData_thisDay$temperature_max, na.rm = TRUE))
  weatherSummary$maxTemperature.max <- c(
    weatherSummary$maxTemperature.max, 
    max(tempData_thisDay$temperature_max, na.rm = TRUE))
  weatherSummary$maxTemperature.min <- c(
    weatherSummary$maxTemperature.min, 
    min(tempData_thisDay$temperature_max, na.rm = TRUE))
  weatherSummary$maxTemperature.error <- c(
    weatherSummary$maxTemperature.error,
    qt(0.975, 
       length(tempData_thisDay$temperature_max) - 1) * 
      sd(tempData_thisDay$temperature_max, na.rm = TRUE) /
      sqrt(length(tempData_thisDay$temperature_max)))
  
  ## daily min
  weatherSummary$minTemperature.mean <- c(
    weatherSummary$minTemperature.mean, 
    mean(tempData_thisDay$temperature_min, na.rm = TRUE))
  weatherSummary$minTemperature.max <- c(
    weatherSummary$minTemperature.max, 
    max(tempData_thisDay$temperature_min, na.rm = TRUE))
  weatherSummary$minTemperature.min <- c(
    weatherSummary$minTemperature.min, 
    min(tempData_thisDay$temperature_min, na.rm = TRUE))
  weatherSummary$minTemperature.error <- c(
    weatherSummary$minTemperature.error,
    qt(0.975, 
       length(tempData_thisDay$temperature_min) - 1) * 
      sd(tempData_thisDay$temperature_min, na.rm = TRUE) /
      sqrt(length(tempData_thisDay$temperature_min)))
  
  ## daily upper and lower deviation
  weatherSummary$temperature.lowerDeviation <- c(
    weatherSummary$temperature.lowerDeviation, 
    mean(tempData_thisDay$temperature - 
           tempData_thisDay$temperature_min))
  weatherSummary$temperature.lowerDeviation.error <- c(
    weatherSummary$temperature.lowerDeviation.error,
    qt(0.975, 
       length(tempData_thisDay$temperature_min) - 1) * 
      sd(tempData_thisDay$temperature - 
           tempData_thisDay$temperature_min, 
         na.rm = TRUE) /
      sqrt(length(tempData_thisDay$temperature_min)))
  weatherSummary$temperature.upperDeviation <- c(
    weatherSummary$temperature.upperDeviation, 
    mean(tempData_thisDay$temperature_max - 
           tempData_thisDay$temperature))
  weatherSummary$temperature.upperDeviation.error <- c(
    weatherSummary$temperature.upperDeviation.error,
    qt(0.975, 
       length(tempData_thisDay$temperature_max) - 1) * 
      sd(tempData_thisDay$temperature_max - 
           tempData_thisDay$temperature, 
         na.rm = TRUE) /
      sqrt(length(tempData_thisDay$temperature_max)))
  
  # precipitation
  weatherSummary$precipitation.mean <- c(
    weatherSummary$precipitation.mean, 
    mean(tempData_thisDay$RAIN, na.rm = TRUE))
  weatherSummary$precipitation.max <- c(
    weatherSummary$precipitation.max, 
    max(tempData_thisDay$RAIN, na.rm = TRUE))
  weatherSummary$precipitation.min <- c(
    weatherSummary$precipitation.min, 
    min(tempData_thisDay$RAIN, na.rm = TRUE))
  weatherSummary$precipitation.error <- c(
    weatherSummary$precipitation.error,
    qt(0.975, 
       length(tempData_thisDay$RAIN) - 1) * 
      sd(tempData_thisDay$RAIN, na.rm = TRUE) /
      sqrt(length(tempData_thisDay$RAIN)))
  
  # ARID
  weatherSummary$ARID.mean <- c(
    weatherSummary$ARID.mean, 
    mean(tempData_thisDay$ARID, na.rm = TRUE))
  weatherSummary$ARID.max <- c(
    weatherSummary$ARID.max, 
    max(tempData_thisDay$ARID, na.rm = TRUE))
  weatherSummary$ARID.min <- c(
    weatherSummary$ARID.min, 
    min(tempData_thisDay$ARID, na.rm = TRUE))
  weatherSummary$ARID.error <- c(
    weatherSummary$ARID.error,
    qt(0.975, 
       length(tempData_thisDay$ARID) - 1) * 
      sd(tempData_thisDay$ARID, na.rm = TRUE) /
      sqrt(length(tempData_thisDay$ARID)))
}
rm(day, tempData_thisDay)

weatherSummary <- data.frame(weatherSummary)
  1. Real-world data
weatherSummary_rakhigarhi <- vector("list", length(weatherSummaryNames))
names(weatherSummary_rakhigarhi) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential

for (day in 1:366)
{
  weatherSummary_rakhigarhi$dayOfYear <- c(weatherSummary_rakhigarhi$dayOfYear, day)
  
  tempData <- weatherData_rakhigarhi[weatherData_rakhigarhi$day == day,]
  
  tempWatbalData <- watbalOutput_rakhigarhi[watbalOutput_rakhigarhi$day == day,]
  
  # solar radiation
  
  weatherSummary_rakhigarhi$solarRadiation.mean <- c(
    weatherSummary_rakhigarhi$solarRadiation.mean, 
    mean(tempData$I, na.rm = TRUE))
  weatherSummary_rakhigarhi$solarRadiation.sd <- c(
    weatherSummary_rakhigarhi$solarRadiation.sd, 
    sd(tempData$I, na.rm = TRUE))
  weatherSummary_rakhigarhi$solarRadiation.max <- c(
    weatherSummary_rakhigarhi$solarRadiation.max, 
    max(tempData$I, na.rm = TRUE))
  weatherSummary_rakhigarhi$solarRadiation.min <- c(
    weatherSummary_rakhigarhi$solarRadiation.min, 
    min(tempData$I, na.rm = TRUE))
  weatherSummary_rakhigarhi$solarRadiation.error <- c(
    weatherSummary_rakhigarhi$solarRadiation.error,
    qt(0.975, length(tempData$I) - 1) * 
      sd(tempData$I, na.rm = TRUE) / 
      sqrt(length(tempData$I)))
  
  # temperature
  
  ## daily mean
  weatherSummary_rakhigarhi$temperature.mean <- c(
    weatherSummary_rakhigarhi$temperature.mean, 
    mean(tempData$T2M, na.rm = TRUE))
  weatherSummary_rakhigarhi$temperature.sd <- c(
    weatherSummary_rakhigarhi$temperature.sd, 
    sd(tempData$T2M, na.rm = TRUE))
  weatherSummary_rakhigarhi$temperature.max <- c(
    weatherSummary_rakhigarhi$temperature.max, 
    max(tempData$T2M, na.rm = TRUE))
  weatherSummary_rakhigarhi$temperature.min <- c(
    weatherSummary_rakhigarhi$temperature.min, 
    min(tempData$T2M, na.rm = TRUE))
  weatherSummary_rakhigarhi$temperature.error <- c(
    weatherSummary_rakhigarhi$temperature.error,
    qt(0.975, length(tempData$T2M) - 1) * 
      sd(tempData$T2M, na.rm = TRUE) / 
      sqrt(length(tempData$T2M)))
  
  ## daily max
  weatherSummary_rakhigarhi$maxTemperature.mean <- c(
    weatherSummary_rakhigarhi$maxTemperature.mean, 
    mean(tempData$Tmax, na.rm = TRUE))
  weatherSummary_rakhigarhi$maxTemperature.max <- c(
    weatherSummary_rakhigarhi$maxTemperature.max, 
    max(tempData$Tmax, na.rm = TRUE))
  weatherSummary_rakhigarhi$maxTemperature.min <- c(
    weatherSummary_rakhigarhi$maxTemperature.min, 
    min(tempData$Tmax, na.rm = TRUE))
  weatherSummary_rakhigarhi$maxTemperature.error <- c(
    weatherSummary_rakhigarhi$maxTemperature.error,
    qt(0.975, length(tempData$Tmax) - 1) * 
      sd(tempData$Tmax, na.rm = TRUE) / 
      sqrt(length(tempData$Tmax)))
  
  ## daily min
  weatherSummary_rakhigarhi$minTemperature.mean <- c(
    weatherSummary_rakhigarhi$minTemperature.mean, 
    mean(tempData$Tmin, na.rm = TRUE))
  weatherSummary_rakhigarhi$minTemperature.max <- c(
    weatherSummary_rakhigarhi$minTemperature.max, 
    max(tempData$Tmin, na.rm = TRUE))
  weatherSummary_rakhigarhi$minTemperature.min <- c(
    weatherSummary_rakhigarhi$minTemperature.min, 
    min(tempData$Tmin, na.rm = TRUE))
  weatherSummary_rakhigarhi$minTemperature.error <- c(
    weatherSummary_rakhigarhi$minTemperature.error,
    qt(0.975, length(tempData$Tmin) - 1) * 
      sd(tempData$Tmin, na.rm = TRUE) / 
      sqrt(length(tempData$Tmin)))
  
  ## daily lower and upper deviation
  weatherSummary_rakhigarhi$temperature.lowerDeviation <- c(
    weatherSummary_rakhigarhi$temperature.lowerDeviation, 
    mean(tempData$T2M - tempData$Tmin)) 
  weatherSummary_rakhigarhi$temperature.lowerDeviation.error <- c(
    weatherSummary_rakhigarhi$temperature.lowerDeviation.error,
    qt(0.975, length(tempData$Tmin) - 1) * 
      sd(tempData$T2M - tempData$Tmin, na.rm = TRUE) / 
      sqrt(length(tempData$Tmin)))
  weatherSummary_rakhigarhi$temperature.upperDeviation <- c(
    weatherSummary_rakhigarhi$temperature.upperDeviation, 
    mean(tempData$Tmax - tempData$T2M))
  weatherSummary_rakhigarhi$temperature.upperDeviation.error <- c(
    weatherSummary_rakhigarhi$temperature.upperDeviation.error,
    qt(0.975, length(tempData$Tmax) - 1) * 
      sd(tempData$Tmax - tempData$T2M, na.rm = TRUE) / 
      sqrt(length(tempData$Tmax)))
  
  # precipitation
  weatherSummary_rakhigarhi$precipitation.mean <- c(
    weatherSummary_rakhigarhi$precipitation.mean, 
    mean(tempData$RAIN, na.rm = TRUE))
  weatherSummary_rakhigarhi$precipitation.max <- c(
    weatherSummary_rakhigarhi$precipitation.max, 
    max(tempData$RAIN, na.rm = TRUE))
  weatherSummary_rakhigarhi$precipitation.min <- c(
    weatherSummary_rakhigarhi$precipitation.min, 
    min(tempData$RAIN, na.rm = TRUE))
  weatherSummary_rakhigarhi$precipitation.error <- c(
    weatherSummary_rakhigarhi$precipitation.error,
    qt(0.975, length(tempData$RAIN) - 1) * 
      sd(tempData$RAIN, na.rm = TRUE) / 
      sqrt(length(tempData$RAIN)))
  
  # ARID
  weatherSummary_rakhigarhi$ARID.mean <- c(
    weatherSummary_rakhigarhi$ARID.mean, 
    mean(tempWatbalData$ARID, na.rm = TRUE))
  weatherSummary_rakhigarhi$ARID.max <- c(
    weatherSummary_rakhigarhi$ARID.max, 
    max(tempWatbalData$ARID, na.rm = TRUE))
  weatherSummary_rakhigarhi$ARID.min <- c(
    weatherSummary_rakhigarhi$ARID.min, 
    min(tempWatbalData$ARID, na.rm = TRUE))
  weatherSummary_rakhigarhi$ARID.error <- c(
    weatherSummary_rakhigarhi$ARID.error,
    qt(0.975, 
       length(tempWatbalData$ARID) - 1) * 
      sd(tempWatbalData$ARID, na.rm = TRUE) /
      sqrt(length(tempWatbalData$ARID)))
}
rm(day, tempData, tempWatbalData)

weatherSummary_rakhigarhi <- data.frame(weatherSummary_rakhigarhi)

Get ranges of weather variables:

roundToMultiple <- function(i, baseOfMultiple, roundFunction = round)
{
  return(match.fun(roundFunction)(i/baseOfMultiple) * baseOfMultiple)
}

rangeSolar = c(
  roundToMultiple(min(c(weatherSummary$solarRadiation.min, weatherSummary_rakhigarhi$solarRadiation.min)), 5, floor),
  roundToMultiple(max(c(weatherSummary$solarRadiation.max, weatherSummary_rakhigarhi$solarRadiation.max)), 5, ceiling))

rangeTemp = c(
  roundToMultiple(min(c(weatherSummary$minTemperature.min, weatherSummary_rakhigarhi$minTemperature.min)), 5, floor),
  roundToMultiple(max(c(weatherSummary$maxTemperature.max, weatherSummary_rakhigarhi$maxTemperature.max)), 5, ceiling))

rangePrecip = c(
  roundToMultiple(min(c(weatherSummary$precipitation.min, weatherSummary_rakhigarhi$precipitation.min)), 5, floor),
  roundToMultiple(max(c(weatherSummary$precipitation.max, weatherSummary_rakhigarhi$precipitation.max)), 5, ceiling))

Set colours for maximum and minimum temperature:

maxTemperatureColour = hsv(7.3/360, 74.6/100, 70/100)

minTemperatureColour = hsv(232/360, 64.6/100, 73/100)

Define auxiliar function for calculating cumulative curve (extracted from weatherModel.R; see documentation folder in Weather model):

getCumulativePrecipitationOfYear <- function(dailyPrecipitationYear)
{
  return(cumsum(dailyPrecipitationYear) / sum(dailyPrecipitationYear))
}

Fig - Weather variables

Create figure:

yearLengthInDays_sim = nlevels(factor(weatherSummary$dayOfYear))
yearLengthInDays_real = nlevels(factor(weatherSummary_rakhigarhi$dayOfYear))

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp1_weatherAndARID.png"
    
    grScale = 2
    fontRescale = 0
    axisTextRescale = 0
    marginTextRescale = 0
    leftPlotMargin = 2
    rightPlotMargin = 4
    cumPrecipMargins <- c(8, 1.5)
    bottomMarginWeight = 2
    cropRangesHeadLength = 0.1
    cropRangesLineWidth = 3
    cropRangesLabelsVerticalSpread = 14
    cropRangesLabelsSize = 0.95
    
    png(plotName, width = grScale * 800, height = grScale * 1000)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp1_weatherAndARID.eps"
    
    grScale = 1.2
    fontRescale = 0.1
    axisTextRescale = -0.1
    marginTextRescale = -0.5
    leftPlotMargin = 1.5
    rightPlotMargin = 3.8
    cumPrecipMargins <- c(6, 1.5)
    bottomMarginWeight = 1
    cropRangesHeadLength = 0.08
    cropRangesLineWidth = 3
    cropRangesLabelsVerticalSpread = 10
    cropRangesLabelsSize = 0.9
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * 8,
                        height = grScale * 10,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  nColumns = 3
  nRowsExceptBottom = 5
  
  layout(rbind(matrix(1:(nColumns * nRowsExceptBottom), 
                      nrow = nRowsExceptBottom,
                      ncol = nColumns, 
                      byrow = FALSE),
               c((nColumns * nRowsExceptBottom) + 1,
                 rep(((nColumns * nRowsExceptBottom) + 2), 2))),
         widths = c(1, 10, 12),
         heights = c(rep(10, nRowsExceptBottom) + c(1, rep(0, nRowsExceptBottom - 2), -2), bottomMarginWeight)
  )
  
  
  
  yLabs <- c(expression(paste("solar radiation (", MJ/m^-2, ")")),
             "temperature (C)", "precipitation (mm)", "ARID", "crop seasons")
  
  par(cex = grScale,
      cex.axis = grScale * (0.8 + axisTextRescale))
  
  # First column: y axis titles
  
  par(mar = c(0, 0, 0, 0.4))
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.5, font = 4, 
       cex = grScale * (0.78 + fontRescale), 
       srt = 90,
       labels = yLabs[1])
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.5, font = 4, 
       cex = grScale * (0.78 + fontRescale), 
       srt = 90,
       labels = yLabs[2])
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.4, font = 4, 
       cex = grScale * (0.78 + fontRescale), 
       srt = 90,
       labels = yLabs[3])
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.6, font = 4, 
       cex = grScale * (0.78 + fontRescale), 
       srt = 90,
       labels = yLabs[4])
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.6, font = 4, 
       cex = grScale * (0.78 + fontRescale), 
       srt = 90,
       labels = yLabs[5])
  
  #==============================================================================
  # Second column: real-world data
  
  # 1. solar radiation
  
  par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
  
  plot(1:yearLengthInDays_real, 
       weatherSummary_rakhigarhi$solarRadiation.mean, 
       axes = FALSE,
       ylim = rangeSolar,
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$solarRadiation.mean + weatherSummary_rakhigarhi$solarRadiation.error),
                rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$solarRadiation.mean - weatherSummary_rakhigarhi$solarRadiation.error),
                rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$solarRadiation.max,
                rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border=NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$solarRadiation.min,
                rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  axis(2, at = seq(rangeSolar[1], rangeSolar[2], 5))
  
  # 2. temperature (daily mean, max, min)
  
  ## daily mean
  plot(1:yearLengthInDays_real, 
       weatherSummary_rakhigarhi$temperature.mean, 
       axes = FALSE,
       ylim = rangeTemp,
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$temperature.mean + weatherSummary_rakhigarhi$temperature.error),
                rev(weatherSummary_rakhigarhi$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$temperature.mean - weatherSummary_rakhigarhi$temperature.error),
                rev(weatherSummary_rakhigarhi$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$temperature.max,
                rev(weatherSummary_rakhigarhi$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$temperature.min,
                rev(weatherSummary_rakhigarhi$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  ## daily max
  lines(1:yearLengthInDays_real, 
        weatherSummary_rakhigarhi$maxTemperature.mean, 
        lwd = grScale, col = maxTemperatureColour)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$maxTemperature.mean + weatherSummary_rakhigarhi$maxTemperature.error),
                rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$maxTemperature.mean - weatherSummary_rakhigarhi$maxTemperature.error),
                rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$maxTemperature.max,
                rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$maxTemperature.min,
                rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
          border = NA)
  
  ## daily min
  lines(1:yearLengthInDays_real, 
        weatherSummary_rakhigarhi$minTemperature.mean, 
        lwd = grScale, col = minTemperatureColour)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$minTemperature.mean + weatherSummary_rakhigarhi$minTemperature.error),
                rev(weatherSummary_rakhigarhi$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$minTemperature.mean - weatherSummary_rakhigarhi$minTemperature.error),
                rev(weatherSummary_rakhigarhi$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$minTemperature.max,
                rev(weatherSummary_rakhigarhi$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$minTemperature.min,
                rev(weatherSummary_rakhigarhi$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  axis(2, at = seq(rangeTemp[1], rangeTemp[2], 5))
  
  # 3. precipitation
  
  # cumulative curve
  par(mar = c(cumPrecipMargins[1], leftPlotMargin, 0.1, 0.1))
  
  plot(c(1, yearLengthInDays_real), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  
  for (year in unique(weatherData_rakhigarhi$year))
  {
    lines(1:nrow(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year,]), 
          getCumulativePrecipitationOfYear(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year, "RAIN"]), 
          lwd = grScale, 
          col = rgb(0, 0, 0, alpha = 0.05))
  }
  rm(year)
  
  # daily values
  par(new = T,
      mar = c(0.1, leftPlotMargin, 0.1, 0.1))
  
  plot(1:yearLengthInDays_real, 
       weatherSummary_rakhigarhi$precipitation.mean, 
       axes = FALSE,
       ylim = rangePrecip * c(1, cumPrecipMargins[2]), 
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$precipitation.mean + weatherSummary_rakhigarhi$precipitation.error),
                rev(weatherSummary_rakhigarhi$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$precipitation.mean - weatherSummary_rakhigarhi$precipitation.error),
                rev(weatherSummary_rakhigarhi$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$precipitation.max,
                rev(weatherSummary_rakhigarhi$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border=NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$precipitation.min,
                rev(weatherSummary_rakhigarhi$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  axis(2, at = seq(rangePrecip[1], rangePrecip[2], 10))
  
  # 4. ARID
  par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
  
  plot(1:yearLengthInDays_real, 
       weatherSummary_rakhigarhi$ARID.mean, 
       axes = FALSE,
       ylim = c(0, 1), 
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$ARID.mean + weatherSummary_rakhigarhi$ARID.error),
                rev(weatherSummary_rakhigarhi$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c((weatherSummary_rakhigarhi$ARID.mean - weatherSummary_rakhigarhi$ARID.error),
                rev(weatherSummary_rakhigarhi$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$ARID.max,
                rev(weatherSummary_rakhigarhi$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border=NA)
  polygon(x = c(1:yearLengthInDays_real, 
                rev(1:yearLengthInDays_real)),
          y = c(weatherSummary_rakhigarhi$ARID.min,
                rev(weatherSummary_rakhigarhi$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  axis(2, at = seq(0, 1, 0.1))
  
  #5. crop growing seasons
  par(mar = c(3, leftPlotMargin, 0.1, 0.1))
  
  plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  for (aCrop in levels(cropTable$crop))
  {
    tempData <- cropTable[cropTable$crop == aCrop,]
    aCropIndex <- match(aCrop, levels(cropTable$crop))
    
    # winter crop
    if (tempData$sugSowingDay > tempData$sugHarvestingDay)
    {
      arrows(x0 = 1, y0 = aCropIndex,
             x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
             angle = 90, code = 3, col = cropColours[aCropIndex],
             length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
      arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
             x1 = 365, y1 = aCropIndex,
             angle = 90, code = 3, col = cropColours[aCropIndex],
             length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
    }
    else
    {
      arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
             x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
             angle = 90, code = 3, col = cropColours[aCropIndex],
             length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
    }
  }
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  # x axis
  axis(1, at = c(0,
                 31, 
                 31+28, 
                 31+28+31, 
                 31+28+31+30, 
                 31+28+31+30+31, 
                 31+28+31+30+31+30,
                 31+28+31+30+31+30+31, 
                 31+28+31+30+31+30+31+31, 
                 31+28+31+30+31+30+31+31+30, 
                 31+28+31+30+31+30+31+31+30+31, 
                 31+28+31+30+31+30+31+31+30+31+30, 
                 31+28+31+30+31+30+31+31+30+31+30+31),
       las = 2
  )
  
  #==============================================================================
  # Third column: simulated data
  
  # 1. solar radiation
  
  par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
  
  plot(1:yearLengthInDays_sim, 
       weatherSummary$solarRadiation.mean, 
       axes = FALSE,
       ylim = rangeSolar,
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$solarRadiation.mean + weatherSummary$solarRadiation.error),
                rev(weatherSummary$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$solarRadiation.mean - weatherSummary$solarRadiation.error),
                rev(weatherSummary$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$solarRadiation.max,
                rev(weatherSummary$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border=NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$solarRadiation.min,
                rev(weatherSummary$solarRadiation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  # 2. temperature (daily mean, max, min)
  
  ## daily mean
  plot(1:yearLengthInDays_sim, 
       weatherSummary$temperature.mean, 
       axes = FALSE,
       ylim = rangeTemp,
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$temperature.mean + weatherSummary$temperature.error),
                rev(weatherSummary$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$temperature.mean - weatherSummary$temperature.error),
                rev(weatherSummary$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$temperature.max,
                rev(weatherSummary$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$temperature.min,
                rev(weatherSummary$temperature.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  ## daily max
  lines(1:yearLengthInDays_sim, 
        weatherSummary$maxTemperature.mean, 
        lwd = grScale, col = maxTemperatureColour)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$maxTemperature.mean + weatherSummary$maxTemperature.error),
                rev(weatherSummary$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$maxTemperature.mean - weatherSummary$maxTemperature.error),
                rev(weatherSummary$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$maxTemperature.max,
                rev(weatherSummary$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$maxTemperature.min,
                rev(weatherSummary$maxTemperature.mean)),
          col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
          border = NA)
  
  ## daily min
  lines(1:yearLengthInDays_sim, 
        weatherSummary$minTemperature.mean, 
        lwd = grScale, col = minTemperatureColour)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$minTemperature.mean + weatherSummary$minTemperature.error),
                rev(weatherSummary$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$minTemperature.mean - weatherSummary$minTemperature.error),
                rev(weatherSummary$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$minTemperature.max,
                rev(weatherSummary$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$minTemperature.min,
                rev(weatherSummary$minTemperature.mean)),
          col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  # 3. precipitation
  
  # cumulative curve
  par(mar = c(cumPrecipMargins[1], 0.1, 0.1, rightPlotMargin))
  
  plot(c(1, yearLengthInDays_sim), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  for (randomSeed in unique(weatherData$randomSeed))
  {
    for (year in unique(weatherData$currentYear))
    {
      lines(1:nrow(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year,]), 
            getCumulativePrecipitationOfYear(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year, "RAIN"]), 
            lwd = grScale, 
            col = rgb(0, 0, 0, alpha = 0.05))
    }
  }
  rm(year, randomSeed)
  
  axis(4, at = seq(0, 1, 0.25))
  mtext("cumulative annual sum", 4, line = 2,
        cex = grScale * (1.5 + marginTextRescale))
  
  # daily values
  par(new = T,
      mar = c(0.1, 0.1, 0.1, rightPlotMargin))
  
  plot(1:yearLengthInDays_sim, 
       weatherSummary$precipitation.mean, 
       axes = FALSE,
       ylim = rangePrecip * c(1, cumPrecipMargins[2]), 
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$precipitation.mean + weatherSummary$precipitation.error),
                rev(weatherSummary$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$precipitation.mean - weatherSummary$precipitation.error),
                rev(weatherSummary$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$precipitation.max,
                rev(weatherSummary$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border=NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$precipitation.min,
                rev(weatherSummary$precipitation.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  # 4. ARID
  par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
  
  plot(1:yearLengthInDays_sim, 
       weatherSummary$ARID.mean, 
       axes = FALSE,
       ylim = c(0, 1), 
       type = "l", lwd = grScale)
  
  ## 95% confidence interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$ARID.mean + weatherSummary$ARID.error),
                rev(weatherSummary$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c((weatherSummary$ARID.mean - weatherSummary$ARID.error),
                rev(weatherSummary$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.5),
          border = NA)
  
  ## min-max interval
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$ARID.max,
                rev(weatherSummary$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border=NA)
  polygon(x = c(1:yearLengthInDays_sim, 
                rev(1:yearLengthInDays_sim)),
          y = c(weatherSummary$ARID.min,
                rev(weatherSummary$ARID.mean)),
          col = rgb(0,0,0, alpha = 0.3),
          border = NA)
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  #5. crop growing seasons
  par(mar = c(3, 0.1, 0.1, rightPlotMargin))
  
  plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  for (aCrop in levels(cropTable$crop))
  {
    tempData <- cropTable[cropTable$crop == aCrop,]
    aCropIndex <- match(aCrop, levels(cropTable$crop))
    
    # winter crop
    if (tempData$sugSowingDay > tempData$sugHarvestingDay)
    {
      arrows(x0 = 1, y0 = aCropIndex,
             x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
             angle = 90, code = 3, col = cropColours[aCropIndex],
             length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
      arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
             x1 = yearLengthInDays_sim, y1 = aCropIndex,
             angle = 90, code = 3, col = cropColours[aCropIndex],
             length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
    }
    else
    {
      arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
             x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
             angle = 90, code = 3, col = cropColours[aCropIndex],
             length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
    }
    
    mtext(aCrop, 
          side = 4, col = cropColours[aCropIndex],
          padj = ((1 - aCropIndex/nrow(cropTable)) - 0.4) * cropRangesLabelsVerticalSpread, 
          las = 2, cex = grScale * cropRangesLabelsSize)
  }
  
  # solstices
  abline(v = 31+28+31+30+31+21, # 21 June (approx.)
         lty = 3, lwd = grScale)
  abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
         lty = 3, lwd = grScale)
  
  # x axis
  axis(1, at = c(0,
                 31, 
                 31+28, 
                 31+28+31, 
                 31+28+31+30, 
                 31+28+31+30+31, 
                 31+28+31+30+31+30,
                 31+28+31+30+31+30+31, 
                 31+28+31+30+31+30+31+31, 
                 31+28+31+30+31+30+31+31+30, 
                 31+28+31+30+31+30+31+31+30+31, 
                 31+28+31+30+31+30+31+31+30+31+30, 
                 31+28+31+30+31+30+31+31+30+31+30+31),
       las = 2
  )
  
  # bottom row: empty and "day of year" or x axis title
  
  par(mar = c(0, 0, 0, 0))
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.45, y = 0.5, font = 4, 
       cex = grScale * (0.8 + fontRescale),
       labels = "day of year")
  
  dev.off()
  
}
rm(imageFormat, grScale, fontRescale, yLabs, nColumns, nRowsExceptBottom, 
   axisTextRescale, marginTextRescale, leftPlotMargin, rightPlotMargin, cumPrecipMargins, bottomMarginWeight,
   tempData, aCrop, aCropIndex, cropRangesHeadLength, cropRangesLineWidth)
knitr::include_graphics(plotName)

Display weather parameters:

knitr::kable(t(yieldData[1, 2:26]), col.names = "default (NW India echo)")
default (NW India echo)
temperature_annualMaxAt2m 37.00
temperature_annualMinAt2m 12.80
temperature_meanDailyFluctuation 2.20
temperature_dailyLowerDeviation 6.80
temperature_dailyUpperDeviation 7.90
solar_annualMax 24.20
solar_annualMin 9.20
solar_meanDailyFluctuation 3.30
CO2_annualMin 245.00
CO2_annualMax 255.00
CO2_meanDailyFluctuation 1.00
precipitation_yearlyMean 489.00
precipitation_yearlySd 142.20
precipitation_dailyCum_nSamples 200.00
precipitation_dailyCum_maxSampleSize 10.00
precipitation_dailyCum_plateauValue_yearlyMean 0.25
precipitation_dailyCum_plateauValue_yearlySd 0.10
precipitation_dailyCum_inflection1_yearlyMean 40.00
precipitation_dailyCum_inflection1_yearlySd 5.00
precipitation_dailyCum_rate1_yearlyMean 0.07
precipitation_dailyCum_rate1_yearlySd 0.02
precipitation_dailyCum_inflection2_yearlyMean 240.00
precipitation_dailyCum_inflection2_yearlySd 20.00
precipitation_dailyCum_rate2_yearlyMean 0.08
precipitation_dailyCum_rate2_yearlySd 0.02

Fig - productivity (yield) per crop type (with default parameter setting)

Summary statistics per crop:

yieldData_summary <- tapply(as.numeric(as.character(yieldData$yield)), yieldData$crop, summary)

yieldData_summary <- data.frame(Reduce(rbind, yieldData_summary), row.names = names(yieldData_summary))
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)

## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)

## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
knitr::kable(yieldData_summary)
Min. X1st.Qu. Median Mean X3rd.Qu. Max.
proso millet 885.6262 1059.3122 1109.8132 1108.5131 1157.0446 1373.1612
pearl millet 300.5779 332.6118 342.6917 342.9984 352.8890 392.0383
rice 0.0000 230.7686 265.1131 261.6839 299.6243 397.6978
barley 337.8926 389.9014 412.2040 414.7265 437.7929 520.7963
wheat 1 280.4535 327.1036 345.8900 347.8794 367.3557 437.3639
wheat 2 296.1396 338.7688 357.1001 359.2555 378.5212 447.0028
for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp1_yieldPerCrop.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 500, height = grScale * 300)
  }
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp1_yieldPerCrop.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 5 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  par(mar = c(6,6,1,1), cex.lab = grScale, cex.axis = 0.8 * grScale)
  
  boxplot(yield ~ factor(crop), data = yieldData,
          ylab = expression(paste("yield (", g/m^2, ")")), xlab = "",
          col = cropColours
          )
  
  dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

ggplot2 option:

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp1_yieldPerCrop_ggplot.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 500, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp1_yieldPerCrop_ggplot.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 5 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  print({
    ggplot(yieldData, aes(y=yield, x=crop, fill=crop)) + 
      geom_violin(position="dodge") +
      #coord_flip() +
      #geom_boxplot() +
      scale_fill_manual(values = cropColours, name="", guide = FALSE) +
      scale_x_discrete() + 
      theme_light(base_size = 11 * grScale) +
      ylab(expression(paste("yield (", g/m^2, ")")))
  })
  dev.off()
}
## Warning: Removed 75 rows containing non-finite values (stat_ydensity).

## Warning: Removed 75 rows containing non-finite values (stat_ydensity).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Fig - mean daily ARID during each crop growing season versus mean ARID of the harvest calendar year

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp1_ARIDVsARID_grow.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 500, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp1_ARIDVsARID_grow.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 5 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  layout(matrix(c(1, 2), nrow = 1, ncol = 2),
         widths = c(10, 3))
  
  par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
  
  boxplot(meanARID_grow ~ factor(crop), data = yieldData,#[yieldData$yield > 0,], # show only non-zero yield 
          ylab = "mean ARID", xlab = "growing seasons",
          #las = 2, 
          col = cropColours)
  
  text(x = 0.5, y = minARID * 1.005, labels = "a", font = 2, cex = grScale)
  
  par(mar = c(5,1,1,1))
  
  xhist <- hist(yieldData$meanARID, breaks = 20, plot = FALSE)
  barplot(xhist$counts, axes = FALSE, space = 0, horiz = TRUE,
          xlab = "current year")
  
  text(x = max(xhist$counts) * 0.9, y = minARID * 1.005, labels = "b", font = 2, cex = grScale)
  
  dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

ggplot2 option:

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp1_ARIDVsARID_grow_ggplot.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 600, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp1_ARIDVsARID_grow_ggplot.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 6 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  tempData <- rbind(yieldData[, c("crop", "meanARID_grow")],
                    cbind(crop = rep("year", nrow(yieldData)), 
                          meanARID_grow = yieldData$meanARID))
  tempData$meanARID_grow <- as.numeric(tempData$meanARID_grow)
  
  print({
    (ggplot(tempData, aes(fill = crop, y = meanARID_grow, x = crop)) + 
       geom_violin(position = "dodge") +
       #geom_boxplot() +
       scale_fill_manual(values = c(cropColours, "grey"), name="", guide = FALSE) +
       scale_x_discrete() + 
       theme_light(base_size = 11 * grScale) +
       theme(axis.title.x = element_blank()) +
       ylab("mean ARID")) +
      ylim(minARID, maxARID)
  })
  
  dev.off()
}
## Warning: Removed 76 rows containing non-finite values (stat_ydensity).

## Warning: Removed 76 rows containing non-finite values (stat_ydensity).
rm(tempData, imageFormat, grScale)
knitr::include_graphics(plotName)

Fig - mean ARID vs productivity (yield)

Plotting absolute values of yield vs ARID:

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp1_yieldVsARID.png"
    
    grScale = 2
    cex.points = 1
    
    png(plotName, width = grScale * 500, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp1_yieldVsARID.eps"
    
    grScale = 2
    cex.points = 0.5
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 5 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
         widths = c(10, 3.5))
  
  par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
  
  plot(c(minARID, maxARID), 
       c(minYield, maxYield),
       xlab = "mean ARID (harvest year)", 
       ylab = expression(paste("yield (", g/m^2, ")")))
  
  for (aCrop in levels(yieldData$crop))
  {
    points(yieldData$meanARID[yieldData$crop == aCrop],
           yieldData$yield[yieldData$crop == aCrop],
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           pch = 20, cex = cex.points
    )
    abline(lm(yield ~ meanARID, data = yieldData[yieldData$crop == aCrop,]),
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           lwd = 2)
  }
  rm(aCrop)
  
  text(x = minARID * 1.005, y = maxYield * 0.95, labels = "a", font = 2, cex = grScale)
  
  plot(c(minARID, maxARID), 
       c(minYield, maxYield),
       xlab = "mean ARID (grow season)",
       ylab = expression(paste("yield (", g/m^2, ")")))
  
  for (aCrop in levels(yieldData$crop))
  {
    points(yieldData$meanARID_grow[yieldData$crop == aCrop],
           yieldData$yield[yieldData$crop == aCrop],
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           pch = 20, cex = cex.points
    )
    abline(lm(yield ~ meanARID_grow, data = yieldData[yieldData$crop == aCrop,]),
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           lwd = 2)
  }
  rm(aCrop)
  
  text(x = minARID * 1.005, y = maxYield * 0.95, labels = "b", font = 2, cex = grScale)
  
  par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  legend(x = 0, y = 1, 
         legend = stringi::stri_c(cropTable$S_water[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"), 
         title = "S_water", 
         fill = cropColours)
  
  dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

ggplot alternative:

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp1_yieldVsARID_ggplot.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 350, height = grScale * 500)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp1_yieldVsARID_ggplot.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 3.5 * grScale,
                          height = 5 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  print({
    
    ggplot(data = yieldData, aes(x = meanARID_grow, y = yield, color = crop, fill = crop)) +
      geom_point(size = 1) +
      geom_smooth(method = 'lm') +
      scale_color_manual(values = cropColours, name="", guide = FALSE) +
      scale_fill_manual(values = cropColours, name="", guide = FALSE) +
      facet_wrap( ~ crop, ncol = 2, scales = 'free') +
      xlab("mean ARID") +
      ylab(expression(paste("yield (", g/m^2, ")"))) +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text = element_text(size = 15),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"))
    
  })
    
  dev.off()
}
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Calculate correlation and linear models fitness:

linearModels.table <- data.frame(
  levels(yieldData$crop),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop))
)

names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")

for (aCrop in levels(yieldData$crop))
{
  linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <- 
    cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop], y = yieldData$yield[yieldData$crop == aCrop], use = "pairwise.complete.obs")
  
  linearModel <- lm(yield ~ meanARID_grow, data = yieldData[yieldData$crop == aCrop,])

  linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
  linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
  linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
  linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
  linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between ARID (x) and Yield (y) per crop")
Pearson’s correlation and linear model between ARID (x) and Yield (y) per crop
crop correlation_Pearson intercept intercept_p speed speed_p adjrsquared
proso millet 0.0201669 1086.7516 0 23.79756 0.5813393 -0.0009296
pearl millet -0.0482022 347.6334 0 -7.44716 0.1872903 0.0009897
rice -0.9094099 693.8981 0 -637.77881 0.0000000 0.8267950
barley -0.8629164 629.6293 0 -273.61862 0.0000000 0.7442715
wheat 1 -0.8470997 531.3303 0 -231.94635 0.0000000 0.7171872
wheat 2 -0.8582940 540.0265 0 -228.55794 0.0000000 0.7363044

Clear temporary objects:

rm(linearModel, linearModels.table, 
   watbalOutput_rakhigarhi, watbalParameters, soil.FC, soil.WP, 
   weatherData, weatherData_rakhigarhi, weatherSummary, weatherSummary_rakhigarhi, weatherSummaryNames,
   xhist, yieldData, yieldData_summary,
   yearLengthInDays_real, yearLengthInDays_sim, 
   maxARID, minARID, maxYield, minYield, rangePrecip, rangeSolar, rangeTemp, 
   maxTemperatureColour, minTemperatureColour, 
   plotName)
## Warning in rm(linearModel, linearModels.table, watbalOutput_rakhigarhi, : object
## 'soil.FC' not found
## Warning in rm(linearModel, linearModels.table, watbalOutput_rakhigarhi, : object
## 'soil.WP' not found

Experiment 2

Experiment 2 perform the same runs of experiment 1, except by varying stochastically two weather model parameters concerning precipitation: precipitation_yearlyMean (average annual sum of precipitation) and precipitation_dailyCum_plateauValue_yearlyMean (average plateau value of the cumulative curve of daily precipitation).

Loading and preparation

Load data:

yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=precipitation-variation_experiment-name=exp2_initRandomSeed=0.csv")

Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):

yieldData$yield[is.na(yieldData$meanARID_grow)] <- NA

Force the same order of crops used by cropTable:

yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))

Calculate the range of ARID, yield, precipitation annual sum (yearly mean) and the “plateau value” (yearly mean) of the cumulative curve of daily precipitation:

minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)

minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)

minPrecipitation_yearlyMean = round(min(yieldData$precipitation_yearlyMean), digits = -2)
maxPrecipitation_yearlyMean = round(max(yieldData$precipitation_yearlyMean), digits = -2)

minPrecipitation_dailyCum_plateauValue_yearlyMean = round(min(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)
maxPrecipitation_dailyCum_plateauValue_yearlyMean = round(max(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)

Display explored values of precipitation_yearlyMean and precipitation_dailyCum_plateauValue_yearlyMean:

layout(matrix(1:2, ncol = 2))
hist(yieldData$precipitation_yearlyMean, main = "precipitation_yearlyMean")
hist(yieldData$precipitation_dailyCum_plateauValue_yearlyMean, main = "precipitation_dailyCum_plateauValue_yearlyMean")

Fig - Effect of precipitation parametres on ARID

Plotting ARID vs precipitation yearly mean and seasonal balance (plateau value):

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp2_ARIDvsPrecipitation.png"
    
    grScale = 2
    cex.points = 1
    
    png(plotName, width = grScale * 500, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp2_ARIDvsPrecipitation.eps"
    
    grScale = 2
    cex.points = 0.5
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 5 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
         widths = c(10, 3.5))
  
  par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
  
  plot(c(minPrecipitation_yearlyMean, maxPrecipitation_yearlyMean),
       c(minARID, maxARID),
       xlab = "precipitation annual sum (mm) (yearly mean)",
       ylab = "mean ARID (growing seasons)")
  
  for (aCrop in levels(yieldData$crop))
  {
    points(yieldData$precipitation_yearlyMean[yieldData$crop == aCrop],
           yieldData$meanARID_grow[yieldData$crop == aCrop],
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           pch = 20, cex = cex.points
    )
    abline(lm(meanARID_grow ~ precipitation_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           lwd = 2)
  }
  rm(aCrop)
  
  text(x = minPrecipitation_yearlyMean * 0.95, y = minARID * 1.2, labels = "a", font = 2, cex = grScale)
  
  plot(c(minPrecipitation_dailyCum_plateauValue_yearlyMean, maxPrecipitation_dailyCum_plateauValue_yearlyMean),
       c(minARID, maxARID),
       xlab = "cumulative daily precipitation plateau value (yearly mean)",
       ylab = "mean ARID (growing seasons)")
  
  for (aCrop in levels(yieldData$crop))
  {
    points(yieldData$precipitation_dailyCum_plateauValue_yearlyMean[yieldData$crop == aCrop],
           yieldData$meanARID_grow[yieldData$crop == aCrop],
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           pch = 20, cex = cex.points
    )
    abline(lm(meanARID_grow ~ precipitation_dailyCum_plateauValue_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           lwd = 2)
  }
  rm(aCrop)
  
  text(x = minPrecipitation_dailyCum_plateauValue_yearlyMean * 0.95, y = minARID * 1.2, labels = "b", font = 2, cex = grScale)
  
  par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  legend(x = 0, y = 1, 
         legend = stringi::stri_c(cropTable$sugSowingDay[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"), 
         title = "sowing day of year", 
         fill = cropColours)
  
  dev.off()
}
rm(imageFormat, grScale, cex.points)
knitr::include_graphics(plotName)

ggplot alternative:

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp2_ARIDvsPrecipitation_ggplot.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 600, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp2_ARIDvsPrecipitation_ggplot.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 5 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  print({
    (ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_yearlyMean, color = crop)) +
       geom_point(size = 1) +
       geom_smooth(method = 'lm') +
       scale_color_discrete_qualitative(guide = FALSE) +
       facet_wrap( ~ crop, ncol = 2) +
       ylab("mean ARID") +
       xlab("precipitation annual sum (mm) (yearly mean)") +
       theme_light(base_size = 11 * grScale) +
       theme(axis.text = element_text(size = 15),
             axis.title = element_text(size = 15),
             strip.background = element_rect(fill = "white"),
             strip.text = element_text(colour = "black"))) +
      (ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_dailyCum_plateauValue_yearlyMean, color = crop)) +
         geom_point(size = 1) +
         geom_smooth(method = 'lm') +
         scale_color_discrete_qualitative(guide = FALSE) +
         facet_wrap( ~ crop, ncol = 2) +
         ylab("") +
         xlab("cumulative daily precipitation plateau value (yearly mean)") +
         theme_light(base_size = 11 * grScale) +
         theme(axis.text = element_text(size = 15),
               axis.title = element_text(size = 15),
               strip.background = element_rect(fill = "white"),
               strip.text = element_text(colour = "black"))) +
      plot_layout(nrow = 1, widths = c(1, 1))
  })
  
  dev.off()
}
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).

## Warning: Removed 75 rows containing missing values (geom_point).
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).

## Warning: Removed 75 rows containing missing values (geom_point).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Calculate correlation and linear models fitness:

linearModels.table <- data.frame(
  levels(yieldData$crop),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop))
)

names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")

for (aCrop in levels(yieldData$crop))
{
  linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <- 
    cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop], 
        y = yieldData$precipitation_yearlyMean[yieldData$crop == aCrop], 
        use = "pairwise.complete.obs")
  
  linearModel <- lm(meanARID_grow ~ precipitation_yearlyMean, data = yieldData[yieldData$crop == aCrop,])

  linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
  linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
  linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
  linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
  linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between precipitation annual sum yearly mean (x) and ARID (y) per crop")
Pearson’s correlation and linear model between precipitation annual sum yearly mean (x) and ARID (y) per crop
crop correlation_Pearson intercept intercept_p speed speed_p adjrsquared
proso millet -0.5760686 0.9960945 0 -0.0001314 0 0.3309618
pearl millet -0.7070320 0.9428775 0 -0.0003964 0 0.4992257
rice -0.7308760 0.9545028 0 -0.0003454 0 0.5335569
barley -0.7133386 0.8883206 0 -0.0004509 0 0.5081726
wheat 1 -0.7138806 0.8916387 0 -0.0004425 0 0.5089473
wheat 2 -0.7138806 0.8916387 0 -0.0004425 0 0.5089473
linearModels.table <- data.frame(
  levels(yieldData$crop),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop))
)

names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")

for (aCrop in levels(yieldData$crop))
{
  linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <- 
    cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop], 
        y = yieldData$precipitation_dailyCum_plateauValue_yearlyMean[yieldData$crop == aCrop], 
        use = "pairwise.complete.obs")
  
  linearModel <- lm(meanARID_grow ~ precipitation_dailyCum_plateauValue_yearlyMean, data = yieldData[yieldData$crop == aCrop,])

  linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
  linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
  linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
  linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
  linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between cumulative daily precipitation plateau value yearly mean (x) and ARID (y) per crop")
Pearson’s correlation and linear model between cumulative daily precipitation plateau value yearly mean (x) and ARID (y) per crop
crop correlation_Pearson intercept intercept_p speed speed_p adjrsquared
proso millet 0.1633846 0.8971123 0 0.0628231 6.9e-06 0.0253933
pearl millet 0.4151101 0.5367639 0 0.3922293 0.0e+00 0.1712098
rice 0.4268577 0.6015239 0 0.3400367 0.0e+00 0.1811142
barley -0.3736394 0.8743004 0 -0.3980747 0.0e+00 0.1384163
wheat 1 -0.3754886 0.8787469 0 -0.3923104 0.0e+00 0.1398036
wheat 2 -0.3754886 0.8787469 0 -0.3923104 0.0e+00 0.1398036

Clear temporary objects:

rm(linearModel, linearModels.table, 
   yieldData,
   maxARID, minARID, maxYield, minYield,
   maxPrecipitation_yearlyMean, minPrecipitation_yearlyMean,
   maxPrecipitation_dailyCum_plateauValue_yearlyMean, minPrecipitation_dailyCum_plateauValue_yearlyMean,
   plotName)

Fig - terrains generated with the Land model

Loading and preparation

Defining terrain random number generator seed (must be available inside “terrains” folder:

listOfTerrainSeeds <- c("0", "35", "56", "72", "92")

Setting dimensions:

figScale = 2
terrainDim = 450 * figScale
topHeight = 50 * figScale
bottomHeight = 150 * figScale
leftWidth = 100 * figScale
internalMargin = 5 * figScale
figWidth = terrainDim * 3 + leftWidth + 2 * internalMargin
figHeight = terrainDim * length(listOfTerrainSeeds) + topHeight + bottomHeight + (length(listOfTerrainSeeds) - 1) * internalMargin

Create separate component images

Create margin header images:

headers <- c("seed",
             "Elevation", 
             "Soil texture",
             "Soil texture type")
namesInFile <- c("seed", "elev", "soilText", "soilTextType")

for (i in 1:length(headers))
{
  #if (i == 1) { width = leftWidth } else { width = terrainDim }
  width = ifelse(i == 1, leftWidth, terrainDim)
  
  png(paste0("FigTerrains_tempFiles/FigTerrains_", namesInFile[i], ".PNG"), 
    width = width, 
    height = topHeight)

  par(mar = rep(0, 4))

  plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, ann = FALSE)
  
  text(0.5, 0.7, cex = figScale * 3.7, adj = c(0.5, 0.5),
       labels = headers[i])
  
  dev.off()
}
rm(i)

for (seed in listOfTerrainSeeds)
{
  png(paste0("FigTerrains_tempFiles/FigTerrains_seed=", seed, ".PNG"), 
    width = leftWidth, 
    height = terrainDim)

  par(mar = rep(0, 4))

  plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
  
  text(0.5, 0.5, cex = figScale * 5, adj = .5,
       labels = seed)
  
  dev.off()
}
rm(seed)

Create legend images:

# code converted from NetLogo (maximum = 255)
elev <- 100/255 + (155/255) * seq(0, 1, length.out = 10)

png("FigTerrains_tempFiles/FigTerrains_legend_elev.PNG", 
    width = terrainDim, 
    height = bottomHeight)

par(mar = c(2,2,1,2), cex = figScale * 2)

barplot(rep(100, 10), 
        col = rgb(elev - (100/255), elev, 0), 
        border = rgb(elev - (100/255), elev, 0),
        axes = FALSE, space = 0) 
mtext("min.", side = 1, line = 0.1, adj = 0, cex = figScale * 2)
mtext("max.", side = 1, line = 0.1, adj = 1, cex = figScale * 2)
#mtext("elevation", side = 1, line = 1.2, adj = .5, cex = figScale * 3)

dev.off()
## png 
##   2
#======================================================================


soilTexture <- c("100% sand", "100% silt", "100% clay")
soilTexture_colors <- c("red", "green", "blue")
xy <- list(c(0.2, 1.05), 
           c(0.2, 0.75), 
           c(0.2, 0.45))

png("FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG", 
    width = terrainDim, 
    height = bottomHeight)

par(mar = rep(0, 4))

plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)

for (i in 1:length(soilTexture))
{
  legend(xy[[i]][1], xy[[i]][2], 
         legend = soilTexture[i], 
         fill = soilTexture_colors[i], 
         bty = "n", adj = 0.1,
         cex = figScale * 2.5, pt.cex = figScale * 5)
}
rm(i)

dev.off()
## png 
##   2
#======================================================================

soilTextureTypes <- c("Sand", "Loamy sand", "Sandy loam", 
                      "Loam", "Silt loam", "Silty clay", 
                      "Clay", "Clay loam", "Sandy clay loam")
soilTextureTypes_colors <- c("#d73229", "#f16a15",  "#9d6e48", 
                             "#eded31", "#59b03c", "#54c4c4", 
                             "#2d8dbe", "#345da9", "#a71b6a")
xy <- list(c(-0.08, 1),   c(0.17,  1), c(0.55,  1), 
           c(-0.08, 0.7), c(0.17, 0.7), c(0.55, 0.7),  
           c(-0.08, 0.4), c(0.17, 0.4), c(0.55, 0.4))

png("FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG", 
    width = terrainDim, 
    height = bottomHeight)

par(mar = rep(0, 4))

plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)

for (i in 1:length(soilTextureTypes))
{
  legend(xy[[i]][1], xy[[i]][2], 
         legend = soilTextureTypes[i], 
         fill = soilTextureTypes_colors[i], 
         bty = "n",  adj = 0.1,
         cex = figScale * 1.8, pt.cex = figScale * 5)
}
rm(i)

dev.off()
## png 
##   2

Create blank image:

png("FigTerrains_tempFiles/FigTerrains_blank.PNG", 
    width = 10, 
    height = 10)
## Warning in png("FigTerrains_tempFiles/FigTerrains_blank.PNG", width = 10, :
## 'width=10, height=10' are unlikely values in pixels
dev.off()
## png 
##   2

Compose final figure

Load input images:

elevTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", 
                       listOfTerrainSeeds ,".PNG")
soilTextTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", 
                           listOfTerrainSeeds ,"_soilTexture.PNG")
soilTextTypeTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", 
                               listOfTerrainSeeds ,"_soilTextureTypes.PNG")

numberOfColumns = 6
numberOfRows = 11

# it is important that these plots were saved as "PNG not "png"
listOfPlots <- c("FigTerrains_tempFiles/FigTerrains_seed.PNG",
                 "FigTerrains_tempFiles/FigTerrains_elev.PNG",
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 "FigTerrains_tempFiles/FigTerrains_soilText.PNG",
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 "FigTerrains_tempFiles/FigTerrains_soilTextType.PNG",
                 # ---
                 paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[1], ".PNG"),
                 elevTerrains[1],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTerrains[1],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTypeTerrains[1],
                 # ---
                 rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
                 # ---
                 paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[2], ".PNG"),
                 elevTerrains[2],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTerrains[2],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTypeTerrains[2],
                 # ---
                 rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
                 # ---
                 paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[3], ".PNG"),
                 elevTerrains[3],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTerrains[3],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTypeTerrains[3],
                 # ---
                 rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
                 # ---
                 paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[4], ".PNG"),
                 elevTerrains[4],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTerrains[4],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTypeTerrains[4],
                 # ---
                 rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
                 # ---
                 paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[5], ".PNG"),
                 elevTerrains[5],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTerrains[5],
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 soilTextTypeTerrains[5],
                 # ---
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 "FigTerrains_tempFiles/FigTerrains_legend_elev.PNG",
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 "FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG",
                 "FigTerrains_tempFiles/FigTerrains_blank.PNG",
                 "FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG"
                 )

grobs <- lapply(lapply(listOfPlots, png::readPNG), grid::rasterGrob)

Create Figure:

widths = 50 * c(leftWidth/figWidth,
                terrainDim/figWidth,
                internalMargin/figWidth,
                terrainDim/figWidth,
                internalMargin/figWidth,
                terrainDim/figWidth) * figScale
heights = (figHeight/figWidth) * 50 * c(topHeight/figHeight, 
                                        terrainDim/figHeight,
                                        internalMargin/figHeight,
                                        terrainDim/figHeight,
                                        internalMargin/figHeight,
                                        terrainDim/figHeight,
                                        internalMargin/figHeight,
                                        terrainDim/figHeight,
                                        internalMargin/figHeight,
                                        terrainDim/figHeight,
                                        bottomHeight/figHeight) * figScale

lay <- matrix(1:(numberOfColumns * numberOfRows), ncol = numberOfColumns, byrow = TRUE)

png("figures/Fig_terrains.png", width = figWidth, height = figHeight)

gridExtra::grid.arrange(grobs = grobs, 
                        layout_matrix = lay,
                        widths = unit(widths, "cm"),
                        heights = unit(heights, "cm"))

dev.off()
## png 
##   2
knitr::include_graphics("figures/Fig_terrains.png")

Show table with parameter values:

terrainsData <- read.csv(
    paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", listOfTerrainSeeds[1] ,".csv"), skip = 8, nrows = 1)

for (i in 2:length(listOfTerrainSeeds))
{
  terrainsData <- rbind(terrainsData, 
                        read.csv(
    paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", listOfTerrainSeeds[i] ,".csv"), skip = 8, nrows = 1))
}
rm(i)
  
terrainsData <- as.data.frame(t(terrainsData[, c(86, 18:39, 91:101, 11:16)]))
names(terrainsData) <- listOfTerrainSeeds
knitr::kable(terrainsData, 
             format = "html",
             col.names = rep("", 5),
             align = c("l", "c", "c", "c"))
randomseed 0 35 56 72 92
elev_algorithm.style “C#” “C#” “C#” “C#” “C#”
elev_featureanglerange 15.866848 1.374065 2.702853 28.400278 25.210475
elev_inversioniterations 1.1360891 1.8251966 1.6684355 0.8325197 0.5830368
elev_noise 3.958625 3.983587 1.630667 2.160736 3.588394
elev_numdepressions 7 10 3 4 6
elev_numprotuberances 4 6 9 4 8
elev_numranges 66 43 38 64 92
elev_numrifts 88 97 13 63 42
elev_rangeaggregation 0.6458941 0.1113466 0.8133660 0.5878475 0.8563670
elev_rangeheight 21.18274020 40.86174048 17.72232293 20.63073219 0.01770265
elev_rangelength 1888 961 680 1279 1659
elev_riftaggregation 0.3834415 0.6789708 0.3916585 0.5006789 0.9385805
elev_riftheight -48.183138 -34.108734 -43.865071 -3.337115 -22.063655
elev_riftlength 3090 959 1589 601 1686
elev_smoothingradius 3.464823 3.464823 3.464823 3.464823 3.464823
elev_smoothstep 1 1 1 1 1
elev_valleyaxisinclination 0.0871293 0.9890636 0.5495040 0.6342591 0.2892803
elev_valleyslope 0.0004043679 0.0037176302 0.0169948110 0.0116752657 0.0080676211
elev_xslope 0.009255966 0.002138160 0.000962116 0.002474697 0.008352354
elev_yslope 0.0007103606 0.0030363731 0.0064754508 0.0040414184 0.0044743707
flow_do.fill.sinks true true true true true
flow_riveraccumulationatstart 42173154 22337132 6731985 13792828 9956658
soil_depthnoise 39.957929 25.789514 35.632002 8.418835 43.506614
soil_formativeerosionrate 2.334470 2.250253 1.322605 1.850611 2.567509
soil_max.clay 95.26008 97.08764 93.12337 11.49762 42.96254
soil_max.sand 88.181043 43.658762 5.469936 84.602521 53.606080
soil_max.silt 68.25092 34.55299 70.62584 76.87028 98.18986
soil_maxdepth 561.0885 149.9677 491.1287 317.9024 464.0493
soil_min.clay 14.33533 76.75279 74.57170 6.94711 35.25457
soil_min.sand 46.147936 39.425832 3.567197 74.532106 17.096166
soil_min.silt 11.82744 33.82258 32.77969 66.05339 77.34002
soil_mindepth 267.5030 105.5965 270.2213 206.9417 229.0143
soil_texturenoise 5.218483 5.335943 9.901911 4.416463 3.734752
ecol_brushfrequencyinflection 13.686815 7.990782 19.641344 7.067255 14.663800
ecol_brushfrequencyrate 0.04661503 0.08916186 0.05192546 0.06610931 0.07624894
ecol_grassfrequencyinflection 2.0733097 0.4220637 4.8174187 3.5072296 0.9123332
ecol_grassfrequencyrate 0.053911121 0.165866239 0.111039067 0.158674553 0.002558042
ecol_woodfrequencyinflection 39.26634 15.82286 35.38110 46.23942 19.50877
ecol_woodfrequencyrate 0.00287898 0.04129579 0.09413197 0.02204882 0.06864615

Clear temporary objects:

rm(grobs, terrainsData,
   elev, soilTexture, soilTexture_colors, xy, 
   elevTerrains, soilTextTerrains, soilTextTypeTerrains,
   soilTextureTypes, soilTextureTypes_colors,
   figScale, terrainDim, topHeight, bottomHeight, leftWidth, internalMargin, figWidth, figHeight, width,
   headers, namesInFile, 
   numberOfColumns, numberOfRows, listOfPlots, 
   widths, heights, lay)

Load the terrain data in the format generated by running experiments 3 and 4 in the I2 model:

terrainData <- list()

for (aTerrainSeed in listOfTerrainSeeds)
{
  #print(paste0("models/output/I2_yield-exp_terrainRandomSeed=", listOfTerrainSeeds[i], ".csv"))
  terrainData[[aTerrainSeed]] <- read.csv(paste0("models/output/I2_yield-exp_terrainRandomSeed=", aTerrainSeed, ".csv"))
}
rm(aTerrainSeed)

Extract map dimensions (using first terrain in listOfTerrainSeeds as reference):

minX = min(terrainData[[1]]$x)
maxX = max(terrainData[[1]]$x)
mapWidth = maxX - minX + 1

minY = min(terrainData[[1]]$y)
maxY = max(terrainData[[1]]$y)
mapHeight = maxY - minY + 1

Calculate the number of land units with river (to be used in terrain classification):

landUnitsWithRiverPerTerrain <- c()

for (aTerrainSeed in listOfTerrainSeeds)
{
  landUnitsWithRiverPerTerrain <- c(landUnitsWithRiverPerTerrain, 
                                    length(which(terrainData[[aTerrainSeed]]$flow_accumulation > mapWidth * mapHeight)))
  # if flow accumulation is greater than the sum of all local contributions, it has a passing river
}
rm(aTerrainSeed)

Experiment 3

Experiment 3 runs the Indus Village’s second integrated model, Land Crop model, using five selected terrains (shown above) and the default parameter configuration aimed at approximating the conditions in Haryana, NW India.

Loading and preparation

Load parameter values:

file_names <- paste0("models/output/", list.files("models/output"))

# load first terrain to initialise yieldData

nameFilter <- grepl(paste0("I2_yield-exp_pars_terrainRandomSeed=", listOfTerrainSeeds[1], 
                           "_type-of-experiment=user-defined_experiment-name=exp3"), file_names)

parsData <- do.call(rbind, lapply(file_names[nameFilter], read.csv))

# load the rest of terrains

for (aTerrainSeed in listOfTerrainSeeds[-1])
{
  nameFilter <- grepl(paste0("I2_yield-exp_pars_terrainRandomSeed=", aTerrainSeed, 
                             "_type-of-experiment=user-defined_experiment-name=exp3"), file_names)
  
  parsData <- rbind(parsData, do.call(rbind, lapply(file_names[nameFilter], read.csv)))
}
rm(aTerrainSeed)

Display all parameters (except cropTable and Land model parameters) default values, as used in experiments 3, 4 and 5:

knitr::kable(t(parsData[1,-c(1, 2)]))
1
temperature_annualMaxAt2m 37
temperature_annualMinAt2m 12.8
temperature_meanDailyFluctuation 2.2
temperature_dailyLowerDeviation 6.8
temperature_dailyUpperDeviation 7.9
solar_annualMax 24.2
solar_annualMin 9.2
solar_meanDailyFluctuation 3.3
precipitation_yearlyMean 489
precipitation_yearlySd 142.2
precipitation_dailyCum_nSamples 200
precipitation_dailyCum_maxSampleSize 10
precipitation_dailyCum_plateauValue_yearlyMean 0.25
precipitation_dailyCum_plateauValue_yearlySd 0.1
precipitation_dailyCum_inflection1_yearlyMean 40
precipitation_dailyCum_inflection1_yearlySd 5
precipitation_dailyCum_rate1_yearlyMean 0.07
precipitation_dailyCum_rate1_yearlySd 0.02
precipitation_dailyCum_inflection2_yearlyMean 240
precipitation_dailyCum_inflection2_yearlySd 20
precipitation_dailyCum_rate2_yearlyMean 0.08
precipitation_dailyCum_rate2_yearlySd 0.02
elev_seaLevelReferenceShift -1000
riverWaterPerFlowAccumulation 1e-04
errorToleranceThreshold 1
crop_selection [‘wheat 1’ ‘wheat 2’ ‘rice’ ‘barley’ ‘pearl millet’ ‘proso millet’]
crop_intensity 50

NOTE: the description of each parameter is writen as comments in the code (I2-integrated-land-crop.nlogo).

Load all experiment 3 output files as a data frame with all terrains:

file_names <- paste0("models/output/", list.files("models/output"))

# load first terrain to initialise yieldData

nameFilter <- grepl(paste0("I2_yield-exp_terrainRandomSeed=", listOfTerrainSeeds[1], 
                           "_type-of-experiment=user-defined_experiment-name=exp3"), file_names)

yieldData <- do.call(rbind, lapply(file_names[nameFilter], read.csv))

# load the rest of terrains

for (aTerrainSeed in listOfTerrainSeeds[-1])
{
  nameFilter <- grepl(paste0("I2_yield-exp_terrainRandomSeed=", aTerrainSeed, 
                             "_type-of-experiment=user-defined_experiment-name=exp3"), file_names)
  
  yieldData <- rbind(yieldData, do.call(rbind, lapply(file_names[nameFilter], read.csv)))
}
rm(aTerrainSeed)

Convert p_crop_yield = 0 to NA when p_soil_meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):

yieldData$p_crop_yield[is.na(yieldData$p_soil_meanARID_grow)] <- NA

Force the same order of crops used by cropTable:

yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))

Assure that random generator seeds are interpreted as factors:

yieldData$terrainRandomSeed <- factor(yieldData$terrainRandomSeed)

yieldData$randomSeed <- factor(yieldData$randomSeed)

Convert production per land unit (total yield) units from g to ton:

yieldData$p_crop_totalYield <- yieldData$p_crop_totalYield / 1E+6

Create factor with river/no river classification of terrains:

hasRiver <- c("no river","river")[factor(landUnitsWithRiverPerTerrain[yieldData$terrainRandomSeed] > 1)]

Fig - productivity (yield) spatial distribution

Build spatial productivity (yield) matrices organised by crop and terrain:

WARNING: this chunk takes a long time to run (i.e. n terrains x 2500 positions x m crops)

meanYield_byPosCropAndTerrain <- list()

for (aTerrainSeed in listOfTerrainSeeds)
{
  meanYield_byPosCropAndTerrain[[aTerrainSeed]] <- list()
  
  # initialise
  
  # total yield
  meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]] <- matrix(rep(NA, mapWidth * mapHeight), nrow = mapHeight, ncol = mapWidth)
  
  for (aCrop in levels(cropTable$crop))
  {
    # yield per crop
    meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]] <- matrix(rep(NA, mapWidth * mapHeight), nrow = mapHeight, ncol = mapWidth)
  }
  
  # fill with values
  
  for (i in minX:maxX)
  {
    for (j in minY:maxY)
    {
      message('\r', paste("terrain", aTerrainSeed, ", position", i, j), appendLF = FALSE)
      
      tempData <- yieldData[yieldData$terrainRandomSeed == aTerrainSeed &        # for this terrain
                              yieldData$x == i & yieldData$y == j,]        # for this position  
      
      tempData2 <- aggregate(tempData$p_crop_yield,                  
                             by = list(randomSeed = tempData$randomSeed,                  # by randomSeed
                                       year = tempData$currentYear),                      # by year
                             FUN = sum)
      
      # mean of the sum of crop yields per year, randomSeed, position and terrainRandomSeed
      meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]][j+1, i+1] <- mean(tempData2$x, na.rm = TRUE)
      
      for (aCrop in levels(cropTable$crop))
      {
        # yield per crop
        meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]][j+1, i+1] <- 
          mean(tempData$p_crop_yield[yieldData$crop == aCrop], na.rm = TRUE)
      }
    }
  }
}
rm(tempData, tempData2,
   i, j, aCrop, aTerrainSeed)

Store spatial matrices in files to ease the workflow (i.e. the chunk above only needs to run if the experiment data has been modified):

for (aTerrainSeed in listOfTerrainSeeds)
{
  # total yield
  write.csv(meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]],
            paste0("FigTerrainYield_tempFiles/exp3_terrainRandomSeed=", aTerrainSeed, "_total.csv"))
  
  for (aCrop in levels(cropTable$crop))
  {
    # yield per crop
    write.csv(meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]],
              paste0("FigTerrainYield_tempFiles/exp3_terrainRandomSeed=", aTerrainSeed, "_", aCrop, ".csv"))
  }
}
rm(aCrop, aTerrainSeed)

Load data from files:

meanYield_byPosCropAndTerrain <- list()

for (aTerrainSeed in listOfTerrainSeeds)
{
  meanYield_byPosCropAndTerrain[[aTerrainSeed]] <- list()
  
  # total yield
  meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]] <- 
    data.matrix(read.csv(paste0("FigTerrainYield_tempFiles/exp3_terrainRandomSeed=", aTerrainSeed, "_total.csv"), row.names = 1))
  
  for (aCrop in levels(cropTable$crop))
  {
    # yield per crop
    meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]] <- 
      data.matrix(read.csv(paste0("FigTerrainYield_tempFiles/exp3_terrainRandomSeed=", aTerrainSeed, "_", aCrop, ".csv"), row.names = 1))
  }
}
rm(aCrop, aTerrainSeed)

Plot mean productivity (mean yield) as grid maps per crop and total (rows) and for each terrain (columns):

mapGridSize = 200 #px
topMargin = mapGridSize / 6
leftMargin = mapGridSize / 2
rightMargin = mapGridSize / 10
bottomMargin = mapGridSize / 2

nColumns = length(meanYield_byPosCropAndTerrain) + 2
nRows = length(meanYield_byPosCropAndTerrain[[1]]) + 2

rowNamesInMargin <- c("seed", gsub(" ", "\n", levels(cropTable$crop)), "total")

leftMarginTextSize = 1.5
topMarginTextSize = 1.2

numBarsInLegend = 10

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp3_yieldSpatialDistribution.png"
    
    grScale = 2
    fontRescale = 0
    
    png(plotName, width = grScale * (mapGridSize * nColumns), height = grScale * (mapGridSize * nRows))
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp3_yieldSpatialDistribution.eps"
    
    grScale = 1.2
    fontRescale = 0.1
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * (mapGridSize * nColumns) / 100,
                        height = grScale * (mapGridSize * nRows) / 100,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  layout(rbind(matrix(1:(nColumns * (nRows - 1)), 
                nrow = nRows - 1, ncol = nColumns, byrow = FALSE),
               rep((nColumns * (nRows - 1)) + 1, nColumns)),
         widths = c(leftMargin, rep(mapGridSize, nColumns - 2), rightMargin),
         heights = c(topMargin, rep(mapGridSize, nRows - 2), bottomMargin))
  
  par(cex = grScale)
  
  # First column: crop names + total
  
  par(mar = c(0, 0, 0, 0.4))
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.5, font = 4, 
       cex = grScale * (topMarginTextSize + fontRescale), 
       labels = rowNamesInMargin[1])
  
  for (aCropIndex in 1:nlevels(cropTable$crop))
  {
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
    text(x = 0.5, y = 0.5, font = 4, 
         cex = grScale * (leftMarginTextSize + fontRescale), 
         #srt = 90,
         labels = rowNamesInMargin[aCropIndex + 1])
  }
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.5, font = 4, 
       cex = grScale * (leftMarginTextSize + fontRescale), 
       #srt = 90,
       labels = rowNamesInMargin[aCropIndex + 2])
  
  #plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  # following columns per each terrain
  
  for (aTerrainSeed in listOfTerrainSeeds)
  {
    # terrain seed
    par(mar = c(0, 0, 0, 0))
    
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
    text(x = 0.5, y = 0.5, font = 4, 
         cex = grScale * (topMarginTextSize + fontRescale), 
         labels = aTerrainSeed)
    
    # map grid for each crop
    par(mar = c(0.1, 0.1, 0.1, 0.1), cex.axis = 0.6 * grScale)
    
    for (aCrop in levels(cropTable$crop))
    {
      # Crop colors can be used instead of the default in image() (i.e. YlOrRd color palette) 
      # colgrp <- findInterval(c(meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]]),
      #                        seq(min(c(meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]])),
      #                            max(c(meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]])),
      #                            length.out = 10))
      # colfunc <- colorRampPalette(c(lighten(cropColours[match(aCrop, levels(cropTable$crop))], amount = 0.9), 
      #                               cropColours[match(aCrop, levels(cropTable$crop))]))
      # collist <- colfunc(length(unique(colgrp))) 
      
      image(t(meanYield_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]]),
            col = hcl.colors(numBarsInLegend, "YlOrRd", rev = TRUE),
            #col = collist, # activate this line to use crop colors
            axes = FALSE, xaxt = 'n', yaxt = 'n', ann = FALSE)
      
      #axis(2, at = seq(0, 1, length.out = mapHeight), labels = minY:maxY)
      #axis(1, at = seq(0, 1, length.out = mapWidth), labels = minX:maxX)
    }
    
    # map grid of the total
    # again, a grey scale can be used intead of the default
    # colgrp <- findInterval(c(meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]]),
    #                        seq(min(c(meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]])),
    #                            max(c(meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]])),
    #                            length.out = 10))
    # colfunc <- colorRampPalette(c(lighten("black", amount = 0.9), "black"))
    # collist <- colfunc(length(unique(colgrp))) 
    
    image(t(meanYield_byPosCropAndTerrain[[aTerrainSeed]][["Total"]]),
          col = hcl.colors(numBarsInLegend, "YlOrRd", rev = TRUE),
          #col = collist, # activate this line to use grey scale colors
          axes = FALSE, xaxt = 'n', yaxt = 'n', ann = FALSE)
  }
  
  # Last column: right margin
  
  for (i in 1:(nRows - 1))
  {
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  }
  
  # bottom margin
  par(mar = c(3,20,1,20), cex = grScale) 
  
  barplot(rep(100, 10), 
          col = hcl.colors(numBarsInLegend, "YlOrRd", rev = TRUE),
          border = hcl.colors(numBarsInLegend, "YlOrRd", rev = TRUE),
          axes = FALSE, space = 0)
  mtext("min.", side = 1, line = 1, adj = 0, cex = grScale * 2)
  mtext("max.", side = 1, line = 1, adj = 1, cex = grScale * 2)  
  
  dev.off()
}
rm(imageFormat, aCropIndex, aCrop, aTerrainSeed, i,
   #colgrp, colfunc, collist,
   mapGridSize, topMargin, leftMargin, rightMargin, bottomMargin, 
   nColumns, nRows, rowNamesInMargin, leftMarginTextSize, topMarginTextSize, numBarsInLegend,
   grScale, fontRescale)
knitr::include_graphics(plotName)

Fig - productivity (yield) per terrain and crop

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = paste0("figures/Fig_exp3_yieldPerTerrainCrop.png")
    
    grScale = 2
    
    png(plotName, width = grScale * 800, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp3_yieldPerTerrainCrop.eps"
    
    grScale = 1.2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * 8,
                        height = grScale * 3,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  print({
    ggplot(yieldData, 
           aes(fill = terrainRandomSeed, y = p_crop_yield, x = crop)) + 
      geom_violin(position="dodge") +
      #geom_boxplot() +
      scale_fill_discrete_diverging(name = "") +
      scale_x_discrete() + 
      facet_wrap( ~ crop, scales = 'free') +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text.x = element_blank(),
            axis.ticks = element_blank(),
            axis.title.x = element_blank(),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"),
            plot.margin = margin(l = 0, unit = "pt")) +
      ylab(expression(paste("yield (", g/m^2, ")")))
  })
  dev.off()
}
## Warning: Removed 375000 rows containing non-finite values (stat_ydensity).

## Warning: Removed 375000 rows containing non-finite values (stat_ydensity).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Fig - productivity (yield) of terrains with and without river

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = paste0("figures/Fig_exp3_yieldPerTerrainCrop_river.png")
    
    grScale = 2
    
    png(plotName, width = grScale * 800, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp3_yieldPerTerrainCrop_river.eps"
    
    grScale = 1.2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * 8,
                        height = grScale * 3,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  print({
    ggplot(cbind(yieldData, extra = hasRiver), 
           aes(fill = extra, y = p_crop_yield, x = extra)) + 
      geom_violin(position="dodge") +
      #geom_boxplot() +
      scale_fill_discrete_diverging(name = "") +
      scale_x_discrete() + 
      facet_wrap( ~ crop, scales = 'free') +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text.x = element_blank(),
            axis.ticks = element_blank(),
            axis.title.x = element_blank(),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"),
            plot.margin = margin(l = 0, unit = "pt")) +
      ylab(expression(paste("yield (", g/m^2, ")")))
  })
  dev.off()
}
## Warning: Removed 375000 rows containing non-finite values (stat_ydensity).

## Warning: Removed 375000 rows containing non-finite values (stat_ydensity).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Clear temporary objects:

rm(file_names, nameFilter, 
   yieldData, hasRiver,
   meanYield_byPosCropAndTerrain, 
   plotName)

Experiment 4

Experiment 4 perform the same runs of experiment 3, except by varying stochastically riverWaterPerFlowAccumulation, representing the water stage increment (mm (height) / m^2 (area)) per unit of flow accumulation at the river’s starting land unit.

Loading and preparation

Load parameter values:

file_names <- paste0("models/output/", list.files("models/output"))

# load first terrain to initialise yieldData

nameFilter <- grepl(paste0("I2_yield-exp_pars_terrainRandomSeed=", listOfTerrainSeeds[1], 
                           "_type-of-experiment=river-variation_experiment-name=exp4"), file_names)

parsData <- do.call(rbind, lapply(file_names[nameFilter], read.csv))

# load the rest of terrains

for (aTerrainSeed in listOfTerrainSeeds[-1])
{
  nameFilter <- grepl(paste0("I2_yield-exp_pars_terrainRandomSeed=", aTerrainSeed, 
                             "_type-of-experiment=river-variation_experiment-name=exp4"), file_names)
  
  parsData <- rbind(parsData, do.call(rbind, lapply(file_names[nameFilter], read.csv)))
}
rm(aTerrainSeed)

Display explored values of riverWaterPerFlowAccumulation:

hist(parsData$riverWaterPerFlowAccumulation, main = "riverWaterPerFlowAccumulation")

knitr::kable(levels(factor(parsData$riverWaterPerFlowAccumulation)))
x
2.02704076590599e-05
8.55452007810128e-05
0.00022977323240978
0.000422851784634783
0.000441634947109688
0.000553325365884843
0.000555289916694427
0.000874695108845879
0.000893931556910296
0.000967359539890516

Load all experiment 4 output files as a data frame with all terrains:

file_names <- paste0("models/output/", list.files("models/output"))

# load first terrain to initialise yieldData

nameFilter <- grepl(paste0("I2_yield-exp_terrainRandomSeed=", listOfTerrainSeeds[1], 
                           "_type-of-experiment=river-variation_experiment-name=exp4"), file_names)

yieldData <- do.call(rbind, lapply(file_names[nameFilter], read.csv))

# load the rest of terrains

for (aTerrainSeed in listOfTerrainSeeds[-1])
{
  nameFilter <- grepl(paste0("I2_yield-exp_terrainRandomSeed=", aTerrainSeed, 
                             "_type-of-experiment=river-variation_experiment-name=exp4"), file_names)
  
  yieldData <- rbind(yieldData, do.call(rbind, lapply(file_names[nameFilter], read.csv)))
}
rm(aTerrainSeed)

Filter out data of land units covered by water surface (i.e. p_ecol_%water = 100 or p_ecol_%crop = 0):

yieldData <- yieldData[yieldData$p_ecol_.water < 100,]

Convert p_crop_yield = 0 to NA when p_soil_meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):

yieldData$p_crop_yield[is.na(yieldData$p_soil_meanARID_grow)] <- NA

Force the same order of crops used by cropTable:

yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))

Assure that random generator seeds are interpreted as factors:

yieldData$terrainRandomSeed <- factor(yieldData$terrainRandomSeed)

yieldData$randomSeed <- factor(yieldData$randomSeed)

Convert production per land unit (total yield) units from g to ton:

yieldData$p_crop_totalYield <- yieldData$p_crop_totalYield / 1E+6

Calculate the range of ARID and productivity (yield):

minARID = round(min(yieldData$p_soil_meanARID_grow, na.rm = TRUE), digits = 2)
maxARID = round(max(yieldData$p_soil_meanARID_grow, na.rm = TRUE), digits = 2)

minYield = round(min(yieldData$p_crop_yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$p_crop_yield, na.rm = TRUE), digits = -1)

minTotalYield = round(min(yieldData$p_crop_totalYield, na.rm = TRUE), digits = 2)
maxTotalYield = round(max(yieldData$p_crop_totalYield, na.rm = TRUE), digits = 2)

minRiverWater = round(min(parsData$riverWaterPerFlowAccumulation), digits = 7)
maxRiverWater = round(max(parsData$riverWaterPerFlowAccumulation), digits = 7)

minPrecipitation_yearlyMean = round(min(parsData$precipitation_yearlyMean), digits = -2)
maxPrecipitation_yearlyMean = round(max(parsData$precipitation_yearlyMean), digits = -2)

minPrecipitation_dailyCum_plateauValue_yearlyMean = round(min(parsData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)
maxPrecipitation_dailyCum_plateauValue_yearlyMean = round(max(parsData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)

Create factor with river/no river classification of terrains:

hasRiver <- c("no river","river")[factor(landUnitsWithRiverPerTerrain[yieldData$terrainRandomSeed] > 1)]

Fig - Effect of river flow on ARID

Plotting ARID vs river flow (only terrains with river):

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp4_ARIDvsRiverFlow.png"
    
    grScale = 2
    cex.points = 1
    
    png(plotName, width = grScale * 500, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp4_ARIDvsRiverFlow.eps"
    
    grScale = 2
    cex.points = 0.5
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 5 * grScale,
                          height = 3 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  layout(matrix(c(1, 2), nrow = 1, ncol = 2, byrow = FALSE),
         widths = c(10, 3.5))
  
  par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
  
  plot(c(minRiverWater, maxRiverWater),
       c(minARID, maxARID),
       xlab = expression(paste("river flow (mm (height) /", m^2, "(area) )")), # ( mm (height) / m^2 (area) )
       ylab = "mean ARID (growing seasons)")
  
  for (aCrop in levels(yieldData$crop))
  {
    tempData <- yieldData[hasRiver == "river" & yieldData$crop == aCrop,]
    
    points(getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                     keyInCalledTable = parsData$randomSeed, 
                                     variableInCalledTable = parsData$riverWaterPerFlowAccumulation),
           tempData$p_soil_meanARID_grow,
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           pch = 20, cex = cex.points
    )
    abline(lm(p_soil_meanARID_grow ~ riverWaterPerFlowAccumulation,
              data = data.frame(cbind(
                p_soil_meanARID_grow = tempData$p_soil_meanARID_grow,
                riverWaterPerFlowAccumulation = getRelatedVariableWithKey(tempData$randomSeed, 
                                                                          parsData$randomSeed, 
                                                                          parsData$riverWaterPerFlowAccumulation)))),
           col = cropColours[match(aCrop, levels(yieldData$crop))],
           lwd = 2)
  }
  
  text(x = minPrecipitation_yearlyMean * 0.95, y = minARID * 1.2, labels = "a", font = 2, cex = grScale)
  
  par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  legend(x = 0, y = 1, 
         legend = stringi::stri_c(cropTable$S_water, " (", levels(yieldData$crop), ")"), 
         title = "S_water (Crop-cultivar)", 
         fill = cropColours)
  
  dev.off()
}
rm(imageFormat, grScale, cex.points, aCrop, tempData)
knitr::include_graphics(plotName)

ggplot alternative:

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp4_ARIDvsRiverFlow_ggplot.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 350, height = grScale * 500)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp4_ARIDvsRiverFlow_ggplot.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 3.5 * grScale,
                          height = 5 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  tempData <- yieldData[hasRiver == "river",]
  
  print({
    ggplot(data = cbind(tempData,
                        riverFlow = getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                                              keyInCalledTable = parsData$randomSeed, 
                                                              variableInCalledTable = parsData$riverWaterPerFlowAccumulation)), 
           aes(x = riverFlow, y = p_soil_meanARID_grow, color = crop, fill = crop)) +
      geom_point(size = 1) +
      geom_smooth(method = 'lm') +
      scale_color_manual(values = cropColours, name="", guide = FALSE) +
      scale_fill_manual(values = cropColours, name="", guide = FALSE) +
      facet_wrap( ~ crop, ncol = 2, scales = 'free') +
      xlab(expression(paste("river flow (mm (height) /", m^2, "(area) )"))) +
      ylab("mean ARID (growing seasons)") +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text = element_text(size = 15),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"))
  })
  dev.off()
}
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 210612 rows containing non-finite values (stat_smooth).
## Warning: Removed 210612 rows containing missing values (geom_point).
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 210612 rows containing non-finite values (stat_smooth).
## Warning: Removed 210612 rows containing missing values (geom_point).
rm(imageFormat, grScale, tempData)
knitr::include_graphics(plotName)

Calculate correlation and linear models fitness:

tempData <- yieldData[hasRiver == "river",]
tempData <- cbind(tempData[,c("crop", "p_soil_meanARID_grow")],
                  riverFlow = getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                                        keyInCalledTable = parsData$randomSeed, 
                                                        variableInCalledTable = parsData$riverWaterPerFlowAccumulation))
linearModels.table <- data.frame(
  levels(yieldData$crop),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop))
)

names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")

for (aCrop in levels(yieldData$crop))
{
  linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <- 
    cor(x = tempData$p_soil_meanARID_grow[tempData$crop == aCrop], 
        y = tempData$riverFlow[tempData$crop == aCrop], 
        use = "pairwise.complete.obs")
  
  linearModel <- lm(p_soil_meanARID_grow ~ riverFlow, data = tempData[tempData$crop == aCrop,])

  linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
  linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
  linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
  linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
  linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop, tempData)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between river flow (x) and ARID (y) per crop")
Pearson’s correlation and linear model between river flow (x) and ARID (y) per crop
crop correlation_Pearson intercept intercept_p speed speed_p adjrsquared
proso millet -0.0440341 0.7828248 0 -49.400547 0.0000000 0.0019361
pearl millet 0.0074461 0.6304191 0 7.870738 0.0000128 0.0000525
rice -0.0011660 0.6589435 0 -1.231364 0.4943451 -0.0000016
barley -0.0425315 0.6815959 0 -49.574053 0.0000000 0.0018053
wheat 1 -0.0428320 0.6850465 0 -49.900830 0.0000000 0.0018309
wheat 2 -0.0428320 0.6850465 0 -49.900830 0.0000000 0.0018309

Fig - Effect of river flow on productivity (yield)

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp4_yieldvsRiverFlow_ggplot.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 350, height = grScale * 500)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp4_yieldvsRiverFlow_ggplot.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 3.5 * grScale,
                          height = 5 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  tempData <- yieldData[hasRiver == "river",]
  
  print({
    ggplot(data = cbind(tempData,
                        riverFlow = getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                                              keyInCalledTable = parsData$randomSeed, 
                                                              variableInCalledTable = parsData$riverWaterPerFlowAccumulation)), 
           aes(x = riverFlow, y = p_crop_yield, color = crop, fill = crop)) +
      geom_point(size = 1) +
      geom_smooth(method = 'lm') +
      scale_color_manual(values = cropColours, name="", guide = FALSE) +
      scale_fill_manual(values = cropColours, name="", guide = FALSE) +
      facet_wrap( ~ crop, ncol = 2, scales = 'free') +
      xlab(expression(paste("river flow (mm (height) /", m^2, "(area) )"))) +
      ylab(expression(paste("yield (", g/m^2, ")"))) +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text = element_text(size = 15),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"))
  })
  dev.off()
}
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 210612 rows containing non-finite values (stat_smooth).
## Warning: Removed 210612 rows containing missing values (geom_point).
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 210612 rows containing non-finite values (stat_smooth).
## Warning: Removed 210612 rows containing missing values (geom_point).
rm(imageFormat, grScale, tempData)
knitr::include_graphics(plotName)

Calculate correlation and linear models fitness:

tempData <- yieldData[hasRiver == "river",]
tempData <- cbind(tempData[,c("crop", "p_crop_yield")],
                  riverFlow = getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                                        keyInCalledTable = parsData$randomSeed, 
                                                        variableInCalledTable = parsData$riverWaterPerFlowAccumulation))
linearModels.table <- data.frame(
  levels(yieldData$crop),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop))
)

names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")

for (aCrop in levels(yieldData$crop))
{
  linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <- 
    cor(x = tempData$p_crop_yield[tempData$crop == aCrop], 
        y = tempData$riverFlow[tempData$crop == aCrop], 
        use = "pairwise.complete.obs")
  
  linearModel <- lm(p_crop_yield ~ riverFlow, data = tempData[tempData$crop == aCrop,])

  linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
  linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
  linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
  linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
  linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(tempData, aCrop)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between river flow (x) and yield (y) per crop")
Pearson’s correlation and linear model between river flow (x) and yield (y) per crop
crop correlation_Pearson intercept intercept_p speed speed_p adjrsquared
proso millet -0.0351944 1113.7076 0 -8784.822 0.0000000 0.0012357
pearl millet -0.0856614 346.1345 0 -5059.249 0.0000000 0.0073350
rice -0.0405137 250.2256 0 -23520.618 0.0000000 0.0016385
barley -0.0074884 437.7860 0 -1925.806 0.0000904 0.0000524
wheat 1 -0.0149060 366.7295 0 -3096.124 0.0000000 0.0002185
wheat 2 -0.0066153 379.9944 0 -1451.674 0.0005431 0.0000401

Effect on total yield:

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp4_totalYieldvsRiverFlow_ggplot.png"
    
    grScale = 2
    
    png(plotName, width = grScale * 350, height = grScale * 500)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp4_totalYieldvsRiverFlow_ggplot.eps"
    
    grScale = 2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::postscript(file = plotName,
                          pointsize = 10,
                          width = 3.5 * grScale,
                          height = 5 * grScale,
                          horizontal = FALSE,
                          paper = "special",
                          onefile = FALSE,
                          family = "sans",
                          colormodel = "cmyk")
  }
  
  tempData <- yieldData[hasRiver == "river",]
  
  print({
    ggplot(data = cbind(tempData,
                        riverFlow = getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                                              keyInCalledTable = parsData$randomSeed, 
                                                              variableInCalledTable = parsData$riverWaterPerFlowAccumulation)), 
           aes(x = riverFlow, y = p_crop_totalYield, color = crop, fill = crop)) +
      geom_point(size = 1) +
      geom_smooth(method = 'lm') +
      scale_color_manual(values = cropColours, name="", guide = FALSE) +
      scale_fill_manual(values = cropColours, name="", guide = FALSE) +
      facet_wrap( ~ crop, ncol = 2, scales = 'free') +
      xlab(expression(paste("river flow (mm (height) /", m^2, "(area) )"))) +
      ylab("total yield (ton) per land unit (ha)") +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text = element_text(size = 15),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"))
  })
  dev.off()
}
## `geom_smooth()` using formula 'y ~ x'
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
## `geom_smooth()` using formula 'y ~ x'
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Calculate correlation and linear models fitness:

tempData <- yieldData[hasRiver == "river",]
tempData <- cbind(tempData[,c("crop", "p_crop_totalYield")],
                  riverFlow = getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                                        keyInCalledTable = parsData$randomSeed, 
                                                        variableInCalledTable = parsData$riverWaterPerFlowAccumulation))
linearModels.table <- data.frame(
  levels(yieldData$crop),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop)),
  numeric(nlevels(yieldData$crop))
)

names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")

for (aCrop in levels(yieldData$crop))
{
  linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <- 
    cor(x = tempData$p_crop_totalYield[tempData$crop == aCrop], 
        y = tempData$riverFlow[tempData$crop == aCrop], 
        use = "pairwise.complete.obs")
  
  linearModel <- lm(p_crop_totalYield ~ riverFlow, data = tempData[tempData$crop == aCrop,])

  linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
  linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
  linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
  linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
  linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop, tempData)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between river flow (x) and yield (y) per crop")
Pearson’s correlation and linear model between river flow (x) and yield (y) per crop
crop correlation_Pearson intercept intercept_p speed speed_p adjrsquared
proso millet -0.2403251 0.9394045 0 -134.18007 0 0.0577534
pearl millet -0.2494387 0.2917571 0 -42.50019 0 0.0622170
rice -0.0407287 0.2085324 0 -19.70730 0 0.0016559
barley -0.0480307 0.2891734 0 -25.42986 0 0.0023040
wheat 1 -0.0509027 0.2422622 0 -22.42551 0 0.0025882
wheat 2 -0.0483261 0.2510331 0 -22.15652 0 0.0023325

Fig - production of land units (total yield) per terrain and crop

In this and the next two figures we can observe the double-edge effect of increased river flow. At one side, more water will be distributed over more land units, stabilising ARID at lower levels throughout the year (i.e. higher yield or productivity). On the other side, however, surface water can accumulate more easily and cover larger extents or the entirety of the area of more land units; particularly those placed in the alluvium, which otherwise are the most productive.

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = paste0("figures/Fig_exp4_totalYieldPerTerrainCrop.png")
    
    grScale = 2
    
    png(plotName, width = grScale * 800, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp4_yieldPerTerrainCrop.eps"
    
    grScale = 1.2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * 8,
                        height = grScale * 3,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  print({
    ggplot(yieldData, 
           aes(fill = terrainRandomSeed, y = p_crop_totalYield, x = crop)) + 
      geom_violin(position="dodge") +
      #geom_boxplot() +
      scale_fill_discrete_diverging(name = "") +
      scale_x_discrete() + 
      facet_wrap( ~ crop, scales = 'free') +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text.x = element_blank(),
            axis.ticks = element_blank(),
            axis.title.x = element_blank(),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"),
            plot.margin = margin(l = 0, unit = "pt")) +
      ylab(expression(paste("yield (", g/m^2, ")")))
  })
  dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

NOTE: remember that terrains 35, 56 and 72 are the ones most affected by changes in river flow.

Fig - production of land units (total yield) of terrains with and without river

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = paste0("figures/Fig_exp4_totalYieldPerTerrainCrop_river.png")
    
    grScale = 2
    
    png(plotName, width = grScale * 800, height = grScale * 300)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp4_totalYieldPerTerrainAndCrop_river.eps"
    
    grScale = 1.2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * 8,
                        height = grScale * 3,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  print({
    ggplot(cbind(yieldData, extra = hasRiver), 
           aes(fill = extra, y = p_crop_totalYield, x = extra)) + 
      geom_violin(position="dodge") +
      #geom_boxplot() +
      scale_fill_discrete_diverging(name = "") +
      scale_x_discrete() + 
      facet_wrap( ~ crop, scales = 'free') +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text.x = element_blank(),
            axis.ticks = element_blank(),
            axis.title.x = element_blank(),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"),
            plot.margin = margin(l = 0, unit = "pt")) +
      ylab("total yield (ton) per land unit (ha)")
  })
  dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Fig - production of land units (total yield) in terrain 35 per crop and explored value of riverWaterPerFlowAccumulation

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = paste0("figures/Fig_exp4_totalYieldPerCropAndRiverflow.png")
    
    grScale = 2
    
    png(plotName, width = grScale * 400, height = grScale * 500)
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp4_totalYieldPerCropAndRiverflow.eps"
    
    grScale = 1.2
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * 4,
                        height = grScale * 5,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  tempData <- yieldData[yieldData$terrainRandomSeed == 35,]
  
  tempData$riverFlow <- getRelatedVariableWithKey(keyInCallingTable =  tempData$randomSeed, 
                                                  keyInCalledTable = parsData$randomSeed, 
                                                  variableInCalledTable = parsData$riverWaterPerFlowAccumulation)
  
  tempData$riverFlow <- factor(formatC(tempData$riverFlow, digits = 2, format = "e"), 
                               levels = formatC(sort(unique(parsData$riverWaterPerFlowAccumulation)), digits = 2, format = "e"))
  
  print({
    ggplot(data = tempData,
           aes(fill = riverFlow, y = p_crop_totalYield, x = riverFlow)) + 
      geom_violin(position = "dodge") +
      #geom_boxplot() +
      scale_fill_discrete_sequential(name = expression(paste("river flow (mm/", m^2, ")")),
                                     guide = guide_legend(reverse = TRUE)) +
      scale_x_discrete() + 
      facet_wrap( ~ crop, scales = 'free', ncol = 1) +
      theme_light(base_size = 11 * grScale) +
      theme(axis.text.x = element_blank(),
            axis.ticks = element_blank(),
            axis.title.x = element_blank(),
            strip.background = element_rect(fill = "white"),
            strip.text = element_text(colour = "black"),
            legend.title = element_text(size = 8 * grScale),
            plot.margin = margin(l = 0, unit = "pt")) +
      ylab("total yield (ton) per land unit (ha)")
  })
  dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)

Clear temporary objects:

rm(file_names, nameFilter,
   linearModel, linearModels.table, 
   yieldData, parsData, hasRiver,
   maxARID, minARID, maxYield, minYield, 
   maxPrecipitation_yearlyMean, minPrecipitation_yearlyMean, 
   maxPrecipitation_dailyCum_plateauValue_yearlyMean, minPrecipitation_dailyCum_plateauValue_yearlyMean,
   maxRiverWater, minRiverWater,
   plotName)

Experiment 5

Experiment 5 perform the same runs of experiment 3, except by varying stochastically the relative frequecies of crops within each land unit. In experiments 3 and 4, frequencies were kept homogeneous and constant, while crop intensity (i.e. the extent of crops in relation to other ecological communities in a land unit) is fixed at 50% throughout all runs. Moreover, in experiments 3 and 4, crop frequency is effectively irrelevant for when analysisng ARID and productivity (yield), as these are measures independent of the relative extent of crops in land units. Frequency is relevant for crop-wise bulk production (total yield or p_crop_totalYield) but the comparisons made remain valid as long as frequencies are homogeneous and constant. The variation of crop frequencies in experiment 5 aims at addressing crop choice factors by evidecing the effect of relative frequencies on bulk production per land unit (i.e. sum of all elements in p_crop_totalYield).

Loading and preparation

Load parameter values:

file_names <- paste0("models/output/", list.files("models/output"))

# load first terrain to initialise yieldData

nameFilter <- grepl(paste0("I2_yield-exp_pars_terrainRandomSeed=", listOfTerrainSeeds[1], 
                           "_type-of-experiment=randomised-crop-frequencies_experiment-name=exp5"), file_names)

parsData <- do.call(rbind, lapply(file_names[nameFilter], read.csv))

# load the rest of terrains

for (aTerrainSeed in listOfTerrainSeeds[-1])
{
  nameFilter <- grepl(paste0("I2_yield-exp_pars_terrainRandomSeed=", aTerrainSeed, 
                             "_type-of-experiment=randomised-crop-frequencies_experiment-name=exp5"), file_names)
  
  parsData <- rbind(parsData, do.call(rbind, lapply(file_names[nameFilter], read.csv)))
}
rm(aTerrainSeed)

Load all experiment 5 output files as a data frame with all terrains:

file_names <- paste0("models/output/", list.files("models/output"))

# load first terrain to initialise yieldData

nameFilter <- grepl(paste0("I2_yield-exp_terrainRandomSeed=", listOfTerrainSeeds[1], 
                           "_type-of-experiment=randomised-crop-frequencies_experiment-name=exp5"), file_names)

yieldData <- do.call(rbind, lapply(file_names[nameFilter], read.csv))

# load the rest of terrains

for (aTerrainSeed in listOfTerrainSeeds[-1])
{
  nameFilter <- grepl(paste0("I2_yield-exp_terrainRandomSeed=", aTerrainSeed, 
                             "_type-of-experiment=randomised-crop-frequencies_experiment-name=exp5"), file_names)
  
  yieldData <- rbind(yieldData, do.call(rbind, lapply(file_names[nameFilter], read.csv)))
}
rm(aTerrainSeed)

Convert p_crop_yield = 0 to NA when p_soil_meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):

yieldData$p_crop_yield[is.na(yieldData$p_soil_meanARID_grow)] <- NA

Force the same order of crops used by cropTable:

yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))

Assure that random generator seeds are interpreted as factors:

yieldData$terrainRandomSeed <- factor(yieldData$terrainRandomSeed)

yieldData$randomSeed <- factor(yieldData$randomSeed)

Convert total yield units (g to ton):

yieldData$p_crop_totalYield <- yieldData$p_crop_totalYield / 1E+6

Calculate the range of yield:

minYield = round(min(yieldData$p_crop_yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$p_crop_yield, na.rm = TRUE), digits = -1)

minTotalYield = round(min(yieldData$p_crop_totalYield), digits = 2)
maxTotalYield = round(max(yieldData$p_crop_totalYield), digits = 2)

Create factor with river/no river classification of terrains:

hasRiver <- c("no river","river")[factor(landUnitsWithRiverPerTerrain[yieldData$terrainRandomSeed] > 1)]

Fig - Crop choice spatial distribution

Plot maps of correlations between total yield of all crops in land unit and the relative frequency of each crop.

Build spatial matrices organised by crop and terrain:

WARNING: this chunk takes a long time to run (i.e. n terrains x 2500 positions x m crops)

corTotalYieldFreq_byPosCropAndTerrain <- list()

for (aTerrainSeed in listOfTerrainSeeds)
{
  corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]] <- list()
  
  # initialise
  
  for (aCrop in levels(cropTable$crop))
  {
    # yield per crop
    corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]] <- matrix(rep(NA, mapWidth * mapHeight), nrow = mapHeight, ncol = mapWidth)
  }
  
  # fill with values
  
  for (i in minX:maxX)
  {
    for (j in minY:maxY)
    {
      message('\r', paste("terrain", aTerrainSeed, ", position", i, j), appendLF = FALSE)
      
      tempData <- yieldData[yieldData$terrainRandomSeed == aTerrainSeed &        # for this terrain
                              yieldData$x == i & yieldData$y == j,]        # for this position  
      
      # sum total yields per year/run
      tempData2 <- aggregate(tempData$p_crop_totalYield,                  
                             by = list(randomSeed = tempData$randomSeed,                  # by randomSeed
                                       year = tempData$currentYear),                      # by year
                             FUN = sum)

      tempData$summedTotalYield <- getRelatedVariableWithTwoKeys(firstKeyInCallingTable = tempData$randomSeed,
                                                                 secondKeyInCallingTable = tempData$currentYear,
                                                                 firstKeyInCalledTable = tempData2$randomSeed,
                                                                 secondKeyInCalledTable = tempData2$year,
                                                                 variableInCalledTable = tempData2$x)
      
      for (aCrop in levels(cropTable$crop))
      {
        tempData3 <- tempData[tempData$crop == aCrop,]
        
        # correlation between frequencies and summed total yield per crop
        corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]][j+1, i+1] <-
          cor(tempData3$p_crop_frequency, tempData3$summedTotalYield)
      }
    }
  }
}
rm(aCrop, i, j, aTerrainSeed,
   tempData, tempData2, tempData3)

Store spatial matrices in files to ease the workflow (i.e. the chunk above only needs to run if the experiment data has been modified):

for (aTerrainSeed in listOfTerrainSeeds)
{
  for (aCrop in levels(cropTable$crop))
  {
    write.csv(corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]],
              paste0("FigTerrainYield_tempFiles/exp5_terrainRandomSeed=", aTerrainSeed, "_", aCrop, ".csv"))
  }
}
rm(aTerrainSeed, aCrop)

Load data from files:

corTotalYieldFreq_byPosCropAndTerrain <- list()

for (aTerrainSeed in listOfTerrainSeeds)
{
  corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]] <- list()
  
  for (aCrop in levels(cropTable$crop))
  {
    # yield per crop
    corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]] <- 
      data.matrix(read.csv(paste0("FigTerrainYield_tempFiles/exp5_terrainRandomSeed=", aTerrainSeed, "_", aCrop, ".csv"), row.names = 1))
  }
}
rm(aTerrainSeed, aCrop)

Calculate range of values crop-wise and terrain-wise:

cropwiseRange <- list()
for (aCrop in levels(cropTable$crop))
{
  tempData <- c()
  for (aTerrainSeed in listOfTerrainSeeds)
  {
    tempData <- c(tempData, as.vector(unlist(corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]])))
  }
  cropwiseRange[[aCrop]] <- c(min(tempData), max(tempData))
}

terrainwiseRange <- list()
for (aTerrainSeed in listOfTerrainSeeds)
{
  tempData <- unlist(corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]])
  
  terrainwiseRange[[aTerrainSeed]] <- c(min(tempData), max(tempData))
}
rm(aCrop, aTerrainSeed, tempData)

Plot mean total yield as grid maps per crop and total (rows) and for each terrain (columns):

mapGridSize = 200 #px
topMargin = mapGridSize / 6
leftMargin = mapGridSize / 2
rightMargin = mapGridSize / 2.5
bottomMargin = mapGridSize / 10

nColumns = length(corTotalYieldFreq_byPosCropAndTerrain) + 2
nRows = length(corTotalYieldFreq_byPosCropAndTerrain[[1]]) + 2

rowNamesInMargin <- c("seed", gsub(" ", "\n", levels(cropTable$crop)))

leftMarginTextSize = 1.5
topMarginTextSize = 1.2

numBarsInLegend = 10

for (imageFormat in listOfImageFormats)
{
  if (imageFormat == "png")
  {
    plotName = "figures/Fig_exp5_corTotalYieldFreqSpatialDistribution.png"
    
    grScale = 2
    fontRescale = 0
    
    png(plotName, width = grScale * (mapGridSize * nColumns), height = grScale * (mapGridSize * nRows))
  }
  
  if (imageFormat == "eps")
  {
    plotName = "figures/Fig_exp5_corTotalYieldFreqSpatialDistribution.eps"
    
    grScale = 1.2
    fontRescale = 0.1
    
    extrafont::loadfonts(device = "postscript")
    grDevices::cairo_ps(filename = plotName,
                        pointsize = 12,
                        width = grScale * (mapGridSize * nColumns) / 100,
                        height = grScale * (mapGridSize * nRows) / 100,
                        onefile = FALSE,
                        fallback_resolution = 600,
                        family = "sans"
    )
  }
  
  layout(rbind(matrix(1:(nColumns * (nRows - 1)), 
                      nrow = nRows - 1, ncol = nColumns, byrow = FALSE),
               ((nColumns * (nRows - 1)) + 1):((nColumns * (nRows - 1)) + nColumns)),
         widths = c(leftMargin, rep(mapGridSize, nColumns - 2), rightMargin),
         heights = c(topMargin, rep(mapGridSize, nRows - 2), bottomMargin))
  
  par(cex = grScale)
  
  # First column: crop names + total
  
  par(mar = c(0, 0, 0, 0.4))
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.5, font = 4, 
       cex = grScale * (topMarginTextSize + fontRescale), 
       labels = rowNamesInMargin[1])
  
  for (aCropIndex in 1:nlevels(cropTable$crop))
  {
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
    text(x = 0.5, y = 0.5, font = 4, 
         cex = grScale * (leftMarginTextSize + fontRescale), 
         #srt = 90,
         labels = rowNamesInMargin[aCropIndex + 1])
  }
  
  # following columns per each terrain
  
  for (aTerrainSeed in listOfTerrainSeeds)
  {
    # terrain seed
    par(mar = c(0, 0, 0, 0))
    
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
    text(x = 0.5, y = 0.5, font = 4, 
         cex = grScale * (topMarginTextSize + fontRescale), 
         labels = aTerrainSeed)
    
    # map grid for each crop
    par(mar = c(0.1, 0.1, 0.1, 0.1), cex.axis = 0.6 * grScale)
    
    for (aCrop in levels(cropTable$crop))
    {
      # Crop colors can be used instead of YlOrRd color palette
      # colgrp <- findInterval(c(corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]]),
      #                        seq(min(c(corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]])),
      #                            max(c(corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]])),
      #                            length.out = 10))
      # colfunc <- colorRampPalette(c(lighten(cropColours[match(aCrop, levels(cropTable$crop))], amount = 0.9), 
      #                               cropColours[match(aCrop, levels(cropTable$crop))]))
      # collist <- colfunc(length(unique(colgrp))) 
      
      image(t(corTotalYieldFreq_byPosCropAndTerrain[[aTerrainSeed]][[aCrop]]),
            zlim = cropwiseRange[[aCrop]], # terrainwiseRange[[aTerrainSeed]],
            col = hcl.colors(numBarsInLegend, "YlOrRd", rev = TRUE),
            #col = collist, # activate this line to use crop colors
            axes = FALSE, xaxt = 'n', yaxt = 'n', ann = FALSE)
      
      #axis(2, at = seq(0, 1, length.out = mapHeight), labels = minY:maxY)
      #axis(1, at = seq(0, 1, length.out = mapWidth), labels = minX:maxX)
    }
  }
  
  # Last column: right margin
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  # scale legend
  par(mar = c(1, 1, 1, 3), cex = grScale)
  
  for (aCrop in levels(cropTable$crop))
  {
    barplot(rep(100, numBarsInLegend), 
            col = hcl.colors(numBarsInLegend, "YlOrRd", rev = TRUE),
            border = hcl.colors(numBarsInLegend, "YlOrRd", rev = TRUE),
            axes = FALSE, space = 0, horiz = TRUE)
    axis(side = 4, at = (0.5):(numBarsInLegend - 0.5),
         labels = round(seq(from = cropwiseRange[[aCrop]][1], to = cropwiseRange[[aCrop]][2], length.out = numBarsInLegend), 1),
         cex.axis = 1, las = 1)
  }
  
  # bottom margin
  
  par(mar = c(0, 0, 0, 0))
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  #par(mar = c(3, 1, 1, 1), cex = grScale) 
  
  for (aTerrainSeed in listOfTerrainSeeds)
  {
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
    # the terrain-wise scale of values can also be shown (remember to apply terrainwise in the zlim argument in image() instead) 
    # barplot(rep(100, numBarsInLegend), 
    #         #names.arg = round(seq(from = min(tempData), to = max(tempData), length.out = numBarsInLegend), 2),
    #         col = heat_hcl(numBarsInLegend, rev = TRUE),
    #         border = heat_hcl(numBarsInLegend, rev = TRUE),
    #         axes = FALSE, space = 0, horiz = FALSE)
    # axis(side = 1, at = 1:numBarsInLegend, 
    #      labels = round(seq(from = terrainwiseRange[[aTerrainSeed]][1], to = terrainwiseRange[[aTerrainSeed]][2], length.out = numBarsInLegend), 2))
  }
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  
  dev.off()
}
rm(imageFormat, aCropIndex, aCrop, aTerrainSeed, i,
   #colgrp, colfunc, collist,
   mapGridSize, topMargin, leftMargin, rightMargin, bottomMargin, 
   nColumns, nRows, rowNamesInMargin, leftMarginTextSize, topMarginTextSize, numBarsInLegend,
   grScale, fontRescale)
## Warning in rm(imageFormat, aCropIndex, aCrop, aTerrainSeed, i, mapGridSize, :
## object 'i' not found
knitr::include_graphics(plotName)

Clear temporary objects:

rm(yieldData, parsData,
   corTotalYieldFreq_byPosCropAndTerrain,
   plotName)